pept-prot-infer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "19f1c6b8-59cb-47b6-9845-be423dd4b230",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "# Summary\n",
    "\n",
    "**Goal: Calculate peptide distribution and protein distribution P(Zj)** from Matt's P(Xi|Zj) fluorosequencing scores (likelihood of a read Xi given peptide Zj)\n",
    "\n",
    "in Juypter Lab: esc, ctrl+a, ctrl+enter to run all cells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa7135a0-23f8-4eef-87d5-d1362653debe",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Abbreviations and so on\n",
    "\n",
    "|Abbreviation|Meaning|Remark|\n",
    "|:-|:-|:-|\n",
    "|P|Likelihood or fraction||\n",
    "|pept|peptide||\n",
    "|prot|protein||\n",
    "|infer|inference||\n",
    "|Z|In pept infer: Pept; in prot infer: protein|E. g. P(Zj) in pept infer is equal to the fraction of each peptide j|\n",
    "|X|In pept infer: Read; in prot infer: protein||\n",
    "|EM|Expectation Maximisation|Algorithm that starts by assuming equal P(Zj) distribution to calculate P(Zj|Xi) and then alternates between updating the two until convergence|\n",
    "|AVG|Average||\n",
    "|STD|Standard deviation||\n",
    "|CI|Confidence interval||"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2277210-a6a0-41e2-b01a-235077e5c406",
   "metadata": {
    "tags": []
   },
   "source": [
    "## What the parts prior to my code do\n",
    "\n",
    "1. **Acquiring of experimental data (reads)**\n",
    "    - 4-5 fluorescence level plots vs Edman degradation cycle, 1 read for each molecule\n",
    "\n",
    "\n",
    "2. **Matts classifier converts reads into P(Xi|Zj)**, the likelihood of each partial sequence given a peptide.\n",
    "    - **Calculating P(Xi) directly impossible/difficult**. You may think that it is simply number of partial sequences/number of all reads, but that is only true if the number of reads is much larger than the number of possible outcomes. It's like rolling a die a couple times, and then asking how likely each number is. You can only do that if you are certain that your number of rolls adequately samples the possibilities. In our case, we do not know this and thus it is instead more consistent to **calculate P(Xi|Zj) instead and to use that to calculate P(Zj), the target variable.**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ae80393-4528-42e8-8236-85d68b550714",
   "metadata": {
    "tags": []
   },
   "source": [
    "## What this code does (well, except for step 2a)\n",
    "\n",
    "**Part 1: Calculation of peptide distribution P(Zj) using Matt's P(Xi|Zj) scores**\n",
    "\n",
    "    1. Imports P(Xi|Zj) table and asks user for settings\n",
    "    2. For each bootstrap run:\n",
    "        1. Create a subtable created from randomly sampling rows in P(Xi|Zj), with replacement (bootstrapping)\n",
    "        2. Calculate P(Zj) via EM using P(Xi|Zj) (subtable)            \n",
    "           - Hen, egg problem. Two missing variables: P(Zj|Xi) and P(Zj). Calculating either requires the other.\n",
    "           - for 1st calculation of P(Zj|Xi) approximate P(Zj) by assuming equal distribution\n",
    "           - alternate between updating P(Zj|Xi) and P(Zj) until convergence\n",
    "\n",
    "    3. Calculate average P(Zj) from all bootstrap runs, and calculate STD, 95% CI**\n",
    "    \n",
    "    \n",
    "**Part 2a: Virtual digest of target proteome and match finding**\n",
    "\n",
    "    1. Virtually digest target proteome into peptides at restriction sites\n",
    "    2. Count number of matches for all peptides inferred from part 1 against digested proteome (save in dyeseq_in_prot_match_count.csv)\n",
    "    3. Calculate P(Xi|Zj) (likelihood of peptide Xi given protein Zj) by dividing the match count of each peptide in a each digested protein by the total number of match counts of ALL peptides against that digested protein\n",
    "    \n",
    "**Part 2b: Calculation of protein distribution using calculated peptide distribution**\n",
    "\n",
    "    1. run bootstrap-EM, like in peptide inferring step\n",
    "        - First calculation of P(Zj|Xi) assumes equal distribution of all P(Zj)\n",
    "        - Uses peptide distribution P(Zj) values from (1) as P(Xi), which simplifies the Bayes-derived equation. However, the principle is the same"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a6677df-1b65-4f1c-bd42-653e4a418d26",
   "metadata": {},
   "source": [
    "## User settings\n",
    "\n",
    "User is asked upon running:\n",
    "\n",
    "|Setting|Meaning|Default|\n",
    "|:-|:-|:-|\n",
    "|EM_convergence_minimum|The lower, the stricter. Checks for difference between prev. and current run|0.0001|\n",
    "|EM_loopcounter_max|Maximum number of EM runs (per bootstrap run)|200|\n",
    "|bootstrap_sampled_fraction|Fraction of subarray sampled for each bootstrap run|0.8|\n",
    "|n_bootstrap_runs|Number of BS runs|200|\n",
    "|CI_percent|Confidence interval|95%|"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2a76c38-258e-4c53-83c2-f00db7b5235d",
   "metadata": {},
   "source": [
    "## Ideas for improvement\n",
    "\n",
    "- Get rid of bootstrap while loop through broadcasting\n",
    "- Get rid of checking for EM convergence every EM run... Instead its probably faster to just e.g. run 199 EM runs, then check for convergence between run 199 and 200. If convergence is still not met, suggest increasing runs to user or just auto-run another 199 runs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0b9c47a-413d-48bb-aba3-421e61fab566",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "080fed80-531c-491a-a919-d15866efa117",
   "metadata": {
    "tags": []
   },
   "source": [
    "## import of packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "b9704e7d-95cf-49e7-816a-9b65de641c1c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np # better arrays than inbuilt arrays\n",
    "import matplotlib.pyplot as plt # to plot stuff\n",
    "\n",
    "import pandas as pd #for DataFrame tables\n",
    "from IPython.display import display #to display dfs more nicely. works similar to head(), but is more flexible in how many columns/rows are shown\n",
    "pd.set_option(\"display.max_rows\", None) #to let display show full df\n",
    "pd.set_option(\"display.max_columns\", None)\n",
    "from sklearn import preprocessing\n",
    "\n",
    "import scipy.stats\n",
    "#from scipy.stats import norm\n",
    "import statistics\n",
    "\n",
    "import math\n",
    "import time #to measure runtime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "307eddb1-7bfa-4c33-aae3-323d2cf41b86",
   "metadata": {
    "tags": []
   },
   "source": [
    "## EM    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "f34fc0bf-52a3-43fb-9269-b40a39654f55",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def update_p_zj_given_xi(p_xi_given_zj, p_zj, p_zj_given_xi):\n",
    "    ## both for PEPTIDE and PROTEIN inference. originally I had thought I could simplify this step for protein inference, but I don't understand why it doesn't work. so the protein inference step uses the same calculcation as the peptide inference step\n",
    "    if inference_mode == \"peptide\" or \"protein\":\n",
    "        numerator = p_xi_given_zj * p_zj # = ARRAY with same shape as P(Xi|Zj)\n",
    "        denominator = np.expand_dims(np.sum(p_xi_given_zj * p_zj, axis=1), axis=-1) # = VECTOR with length of p Xi. expand_dims is necessary because else denominator has the shape (number of reads Xi,) instead of (number of reads Xi, 1). The entries of the vector are identical. But it does not have the same matrix \"rank\", which is why broadcasting fails without that \n",
    "        p_zj_given_xi = np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator!=0) # if denominator = 0, return 0. This happens sometimes when the likelihood for Zj is equal to 0. This could alternatively be avoided by cropping the rows where the input P(Xi|Zj) scores equal to 0\n",
    "        p_zj_given_xi = np.transpose(p_zj_given_xi)\n",
    "       \n",
    "    ## UNUSED. PLACEHOLDER for PROTEIN inference (denominator is has no loop). See inference_mode comment above for explanation why its unused.\n",
    "    elif inference_mode == \"PLACEHOLDER protein\":\n",
    "        p_zj_given_xi = (p_xi_given_zj * p_zj)/np.expand_dims(p_xi, axis=-1)\n",
    "        p_zj_given_xi = np.transpose(p_zj_given_xi)\n",
    "        \n",
    "    else:\n",
    "        print(\"ERROR: inference_mode was not set to peptide or protein\")\n",
    "        \n",
    "    return p_zj_given_xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "d0dfeb6f-3479-4489-8502-8ec33cab074e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def update_p_zj(p_zj, p_zj_given_xi):  \n",
    "    if inference_mode == \"peptide\":\n",
    "        p_zj = p_zj_given_xi.sum(axis=1)#/p_xi_given_zj.shape[1]\n",
    "\n",
    "        p_zj_sum = np.sum(p_zj, axis = 0)\n",
    "        p_zj = p_zj/p_zj_sum #so that p_zj is a fraction (i. e. a number between 0 and 1) \n",
    "\n",
    "    elif inference_mode == \"protein\":\n",
    "        p_zj.fill(0) # reset all values of pzj to 0\n",
    "        \n",
    "        p_zj = np.sum(p_zj_given_xi * p_xi, axis = 1)\n",
    "        p_zj_sum = np.sum(p_zj, axis = 0)\n",
    "        p_zj = p_zj/p_zj_sum #so that p_zj is a fraction (i. e. a number between 0 and 1) \n",
    "    \n",
    "    else:\n",
    "        print(\"ERROR: inference_mode was not set to peptide or protein\")\n",
    "    \n",
    "    return p_zj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "3c5b5972-b626-454b-8422-c5be7e306557",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def EM_convergence_checker(p_zj_old, p_zj, EM_convergence_minimum):\n",
    "    difference_abs = abs(np.sum(p_zj_old[round(0.5*len(p_zj_old))::], axis=0, dtype=float) - np.sum(p_zj[round(0.5*len(p_zj_old))::], axis=0, dtype=float))\n",
    "    if difference_abs > EM_convergence_minimum:\n",
    "        return False\n",
    "    \n",
    "    else:\n",
    "        return True   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "79e08316-e78a-4de3-b949-3ba2077fd813",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def EM(p_xi_given_zj):\n",
    "    global xi_len\n",
    "    global zj_len\n",
    "    xi_len = p_xi_given_zj.shape[0] # length of x (number of reads in pept infer, or peptides in prot infer)\n",
    "    zj_len = p_xi_given_zj.shape[1] # length of z (number of peptides in pept infer, or proteins in prot infer)\n",
    "    \n",
    "    p_zj_given_xi = np.full((zj_len, xi_len), 0, dtype=float) #Initialisation based on array size of p_xi_given_zj -- same size, but transposed\n",
    "    \n",
    "    p_zj_initial = 1/zj_len #initial approximation: all zj equally likely, to jumpstart first iteration\n",
    "\n",
    "    global p_zj\n",
    "    p_zj = np.full(zj_len, p_zj_initial)\n",
    "    p_zj_old = np.full(zj_len, 0)\n",
    "    \n",
    "    global EM_loopcounter \n",
    "    EM_loopcounter = 0\n",
    "    EM_convergence_checker_result = False\n",
    "        \n",
    "    while (EM_convergence_checker_result == False) and EM_loopcounter < EM_loopcounter_max:\n",
    "        p_zj_given_xi = update_p_zj_given_xi(p_xi_given_zj, p_zj, p_zj_given_xi)        \n",
    "        p_zj_old = np.copy(p_zj, order='K', subok=False)\n",
    "        p_zj = update_p_zj(p_zj, p_zj_given_xi)\n",
    "        EM_convergence_checker_result = EM_convergence_checker(p_zj_old, p_zj, EM_convergence_minimum)\n",
    "\n",
    "        EM_loopcounter = EM_loopcounter + 1\n",
    "    return p_zj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5084a15-17a8-498a-a95e-024edc219205",
   "metadata": {
    "tags": []
   },
   "source": [
    "## bootstrapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "017f6966-6cfe-47b9-b11c-46cbe649985e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def create_subarray_of_p_xi_given_zj(p_xi_given_zj):\n",
    "    df_p_xi_given_zj = pd.DataFrame(p_xi_given_zj)\n",
    "    df_p_xi_given_zj_sample = df_p_xi_given_zj.sample(frac=bootstrap_sampled_fraction, axis='rows', replace=True) # filters for a random partial dataset\n",
    "    # display(pd.DataFrame(df_p_xi_given_zj_sample))\n",
    "    p_xi_given_zj_subarray = df_p_xi_given_zj_sample.to_numpy()\n",
    "    \n",
    "    # print(\"p_xi_given_zj_subarray\")\n",
    "    # display(pd.DataFrame(p_xi_given_zj_subarray))\n",
    "    \n",
    "    return p_xi_given_zj_subarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "4edffbe3-83fa-4e42-bf2f-17beafeb014a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def bootstrap_EM():\n",
    "    global p_zj_bootstrap_results_fraction\n",
    "    global p_zj_peptide_copy\n",
    "    global p_xi\n",
    "    i = 0\n",
    "    \n",
    "    p_zj_bootstrap_results_absolute = np.full((p_xi_given_zj.shape[1]), 0, dtype=float)\n",
    "      \n",
    "    if bootstrap_sampled_fraction == -1: # no bootstrapping, i. e. use full dataset for EM. Since the EM is deterministic, there is no point in running multiple bootstrap runs\n",
    "        if inference_mode == \"peptide\": \n",
    "            p_xi_given_zj_subarray = p_xi_given_zj\n",
    "            p_zj_bootstrap_results_absolute = EM(p_xi_given_zj_subarray)\n",
    "\n",
    "            p_zj_bootstrap_results_fraction = p_zj_bootstrap_results_absolute/np.sum(p_zj_bootstrap_results_absolute, axis = 0)\n",
    "            print(\"Peptide Zj values from every bootstrap run (columns: peptides, rows: bootstrap run, displayed as fractions: \")\n",
    "    \n",
    "        elif inference_mode == \"protein\": #only difference: uses p(zj) values from peptide inferrence as the p_xi values. Also, update_p_zj() called by EM() uses a different formula than in pep infer\n",
    "            p_xi = np.copy(p_zj_bootstrap_results_fraction_avg) # input is just a vector, because no bootstrapping.\n",
    "\n",
    "            p_xi_given_zj_subarray = p_xi_given_zj\n",
    "            p_zj_bootstrap_results_absolute = EM(p_xi_given_zj_subarray)\n",
    "\n",
    "            p_zj_bootstrap_results_fraction = p_zj_bootstrap_results_absolute/np.sum(p_zj_bootstrap_results_absolute, axis = 0)\n",
    "            print(\"PRE-CORRECTION: Protein Zj values from every bootstrap run (columns: protein, rows: bootstrap run, displayed as fractions: \")\n",
    "    \n",
    "    elif bootstrap_sampled_fraction != -1: # with bootstrapping, i. e. a subarray of the dataset is created for each bootstrap run\n",
    "        if inference_mode == \"peptide\":\n",
    "            while i < n_bootstrap_runs:\n",
    "                p_xi_given_zj_subarray = create_subarray_of_p_xi_given_zj(p_xi_given_zj) \n",
    "                if i == 0:\n",
    "                    p_zj_bootstrap_results_absolute = EM(p_xi_given_zj_subarray)\n",
    "                else:\n",
    "                    p_zj_bootstrap_results_absolute = np.vstack((p_zj_bootstrap_results_absolute, EM(p_xi_given_zj_subarray)))\n",
    "\n",
    "                print(\"Bootstrap run #\", i, \". EM loops: \", EM_loopcounter, sep=\"\")\n",
    "                if EM_loopcounter == EM_loopcounter_max:\n",
    "                    print(\"WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\")\n",
    "\n",
    "                i = i + 1\n",
    "            p_zj_bootstrap_results_fraction = p_zj_bootstrap_results_absolute/np.sum(p_zj_bootstrap_results_absolute, axis = 1)[0]    \n",
    "\n",
    "        elif inference_mode == \"protein\":\n",
    "            p_zj_peptide_copy = np.copy(p_zj_bootstrap_results_fraction) # input: an array\n",
    "\n",
    "            while i < n_bootstrap_runs:\n",
    "                p_xi_given_zj_subarray = p_xi_given_zj\n",
    "                p_xi = p_zj_peptide_copy[i]\n",
    "                \n",
    "                if i == 0:\n",
    "                    p_zj_bootstrap_results_absolute = EM(p_xi_given_zj_subarray)\n",
    "                else:\n",
    "                    p_zj_bootstrap_results_absolute = np.vstack((p_zj_bootstrap_results_absolute, EM(p_xi_given_zj_subarray)))\n",
    "\n",
    "                print(\"Bootstrap run #\", i, \". EM loops: \", EM_loopcounter, sep=\"\")\n",
    "                if EM_loopcounter == EM_loopcounter_max:\n",
    "                    print(\"WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\")\n",
    "\n",
    "                i = i + 1\n",
    "                \n",
    "            p_zj_bootstrap_results_fraction = p_zj_bootstrap_results_absolute/np.sum(p_zj_bootstrap_results_absolute, axis = 1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "9df928ed-adcb-402d-9ec5-97ffd317a4d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def bootstrap_EM_analytics_AVG():\n",
    "    global p_zj_bootstrap_results_fraction_avg\n",
    "    \n",
    "    if bootstrap_sampled_fraction == -1:\n",
    "        p_zj_bootstrap_results_fraction_avg = np.copy(p_zj_bootstrap_results_fraction,order='K')\n",
    "    elif bootstrap_sampled_fraction != -1:\n",
    "        p_zj_bootstrap_results_fraction_avg = np.sum(p_zj_bootstrap_results_fraction, axis = 0)/p_zj_bootstrap_results_fraction.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "ec5ce1eb-7e35-4cf1-b20a-b632273a4487",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def bootstrap_EM_analytics_CI():\n",
    "    outlier_on_each_side = (100-CI_percent)/(2*100) # i.e. 2.5% missing on each side if 95% CI\n",
    "    \n",
    "    lowindex = round(outlier_on_each_side*n_bootstrap_runs)\n",
    "    highindex = round((1-outlier_on_each_side)*n_bootstrap_runs)\n",
    "    \n",
    "    bootstrap_CI_max = np.sort(p_zj_bootstrap_results_fraction, axis=0, kind=\"quicksort\", order=None)[lowindex:highindex:][-1] #fetches last row from sorted table, ie the highest values. axis=0: sorting along each column\n",
    "    bootstrap_CI_min = np.sort(p_zj_bootstrap_results_fraction, axis=0, kind=\"quicksort\", order=None)[lowindex:highindex:][0] #same but for first row, highest values\n",
    "    bootstrap_CI_minmax = np.stack((bootstrap_CI_min, bootstrap_CI_max))\n",
    "    \n",
    "    ## show full sorted table\n",
    "    # global p_zj_bootstrap_results_fraction_CI\n",
    "    # print(\"The\", CI_percent, \"% confidence interval of all zj bootstrapping values sorted, displayed as fractions: \")    \n",
    "    # p_zj_bootstrap_results_fraction_CI = np.sort(p_zj_bootstrap_results_fraction, axis=0, kind=\"quicksort\", order=None)[lowindex:highindex:] #axis=sorting along each column\n",
    "    # display(pd.DataFrame(p_zj_bootstrap_results_fraction_CI))\n",
    "    \n",
    "    return bootstrap_CI_minmax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0d7e608-f5e3-4c2e-af04-04736ece38fe",
   "metadata": {
    "tags": []
   },
   "source": [
    "## output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "d65e8230-c958-4e74-98df-de37b323b21b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def output():\n",
    "    global p_zj_bootstrap_results_fraction_avg\n",
    "    bootstrap_EM_analytics_AVG()\n",
    "    bootstrap_CI_minmax = bootstrap_EM_analytics_CI()\n",
    "\n",
    "    ### STDs (alternative to CI)\n",
    "    p_zj_bootstrap_results_fraction_std = np.std(p_zj_bootstrap_results_fraction, axis = 0)\n",
    "    \n",
    "    if inference_mode == \"protein\":\n",
    "        corrfactor = dyeseq_in_prot_match_count*np.sum(p_zj_bootstrap_results_fraction_avg/dyeseq_in_prot_match_count)\n",
    "\n",
    "        p_zj_bootstrap_results_fraction_avg = p_zj_bootstrap_results_fraction_avg/corrfactor        \n",
    "        if bootstrap_sampled_fraction != -1:\n",
    "            p_zj_bootstrap_results_fraction_std = p_zj_bootstrap_results_fraction_std/corrfactor\n",
    "            bootstrap_CI_minmax = bootstrap_CI_minmax/corrfactor\n",
    "            \n",
    "    # Combines AVG zj values (from all bootstrap runs), the STD, and the two bounds for the user-chosen confidence interval in one  \n",
    "    if bootstrap_sampled_fraction == -1: #ie if no bootstrapping. Fills -CI and +CI columns with \"N/A\" values\n",
    "        p_zj_avg_plus_CI = np.stack((p_zj_bootstrap_results_fraction_avg*100, np.full((len(p_zj_bootstrap_results_fraction_avg)), fill_value=\"N/A\", dtype=None), np.full((len(p_zj_bootstrap_results_fraction_avg)), fill_value=\"N/A\", dtype=None), np.full((len(p_zj_bootstrap_results_fraction_avg)), fill_value=\"N/A\", dtype=None)), axis=1)    \n",
    "        display(pd.DataFrame(p_zj_avg_plus_CI, columns = [\"AVG [%]\", \"±STD [%]\", \"-CI [%]\", \"+CI [%]\"]).astype({\"AVG [%]\":float}).round(4))\n",
    "    else:\n",
    "        p_zj_avg_plus_CI = np.stack((p_zj_bootstrap_results_fraction_avg, p_zj_bootstrap_results_fraction_std, bootstrap_CI_minmax[0], bootstrap_CI_minmax[1]), axis = 1)\n",
    "        display(pd.DataFrame(p_zj_avg_plus_CI, columns = [\"AVG [%]\", \"±STD [%]\", \"-CI [%]\", \"+CI [%]\"]).round(4)*100)\n",
    "\n",
    "    ### print user-settings\n",
    "    print('\\033[1m' + 'Settings of this run' + '\\033[0m')\n",
    "    print(\"Inference mode:\", inference_mode)\n",
    "    print(\"EM_convergence_minimum:\", EM_convergence_minimum)\n",
    "    print(\"EM_loopcounter_max:\", EM_loopcounter_max)\n",
    "    print(\"bootstrap_sampled_fraction:\", bootstrap_sampled_fraction)\n",
    "    print(\"n_bootstrap_runs:\", n_bootstrap_runs)\n",
    "    print(\"CI_percent:\", CI_percent, \"\\n\")\n",
    "\n",
    "    ### print other\n",
    "    print(\"Runtime: %s seconds ---\" % round((time.time() - start_time),2)) #run time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d8adb5e-3057-424e-866e-637ebccd86ee",
   "metadata": {
    "tags": []
   },
   "source": [
    "# call functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "c6c157ec-dcbd-457a-a448-bde516d2728d",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdin",
     "output_type": "stream",
     "text": [
      "EM_convergence_minimum? If nothing is entered, it is set to 0.0001. \n",
      "Maximum number of EM runs (per bootstrap run)? If nothing is entered, it is set to 200. 10\n",
      "Fraction of subarray sampled for each bootstrap run? If nothing is entered, it is set to 0.8. If -1 is entered, bootstrapping is turned off. Note: Bootstrapping is always turned OFF for the protein inference part. \n",
      "Number of bootstrap runs? If nothing is entered, it is set to 200. 10\n",
      "Condidence interval? If nothing is entered, it is set to 95. \n"
     ]
    }
   ],
   "source": [
    "# esc, ctrl+a, ctrl+enter to run all cells\n",
    "\n",
    "EM_convergence_minimum = float(input(\"EM_convergence_minimum? If nothing is entered, it is set to 0.0001.\") or \"0.0001\")\n",
    "EM_loopcounter_max = int(input(\"Maximum number of EM runs (per bootstrap run)? If nothing is entered, it is set to 200.\") or \"200\") \n",
    "bootstrap_sampled_fraction = float(input(\"Fraction of subarray sampled for each bootstrap run? If nothing is entered, it is set to 0.8. If -1 is entered, bootstrapping is turned off. Note: Bootstrapping is always turned OFF for the protein inference part.\") or \"0.8\")                   \n",
    "n_bootstrap_runs = int(input(\"Number of bootstrap runs? If nothing is entered, it is set to 200.\") or \"200\") \n",
    "CI_percent = float(input(\"Condidence interval? If nothing is entered, it is set to 95.\") or \"95\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "6fdc4f76-547b-4355-8966-3d17d772a815",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bootstrap run #0. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #1. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #2. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #3. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #4. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #5. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #6. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #7. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #8. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #9. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n"
     ]
    },
    {
     "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>AVG [%]</th>\n",
       "      <th>±STD [%]</th>\n",
       "      <th>-CI [%]</th>\n",
       "      <th>+CI [%]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.11</td>\n",
       "      <td>0.12</td>\n",
       "      <td>2.87</td>\n",
       "      <td>3.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.66</td>\n",
       "      <td>0.09</td>\n",
       "      <td>4.53</td>\n",
       "      <td>4.79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.58</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.50</td>\n",
       "      <td>0.65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.55</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.47</td>\n",
       "      <td>0.71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.54</td>\n",
       "      <td>0.13</td>\n",
       "      <td>2.25</td>\n",
       "      <td>2.70</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>1.19</td>\n",
       "      <td>0.06</td>\n",
       "      <td>1.10</td>\n",
       "      <td>1.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>1.57</td>\n",
       "      <td>0.12</td>\n",
       "      <td>1.40</td>\n",
       "      <td>1.81</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4.62</td>\n",
       "      <td>0.11</td>\n",
       "      <td>4.41</td>\n",
       "      <td>4.76</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>1.24</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.98</td>\n",
       "      <td>1.37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.78</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.70</td>\n",
       "      <td>0.92</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.41</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.36</td>\n",
       "      <td>0.47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0.59</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.52</td>\n",
       "      <td>0.63</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.15</td>\n",
       "      <td>0.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>2.23</td>\n",
       "      <td>0.11</td>\n",
       "      <td>2.06</td>\n",
       "      <td>2.43</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0.93</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.84</td>\n",
       "      <td>1.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>4.62</td>\n",
       "      <td>0.17</td>\n",
       "      <td>4.29</td>\n",
       "      <td>5.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.73</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.69</td>\n",
       "      <td>0.84</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>0.40</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.34</td>\n",
       "      <td>0.45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.68</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.61</td>\n",
       "      <td>0.78</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>0.53</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.45</td>\n",
       "      <td>0.59</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>0.09</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>0.55</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.41</td>\n",
       "      <td>0.66</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>2.60</td>\n",
       "      <td>0.12</td>\n",
       "      <td>2.38</td>\n",
       "      <td>2.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>3.23</td>\n",
       "      <td>0.11</td>\n",
       "      <td>3.07</td>\n",
       "      <td>3.40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>0.08</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>1.79</td>\n",
       "      <td>0.06</td>\n",
       "      <td>1.67</td>\n",
       "      <td>1.88</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>0.27</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>0.62</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.54</td>\n",
       "      <td>0.77</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>0.26</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>0.32</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.26</td>\n",
       "      <td>0.37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>3.32</td>\n",
       "      <td>0.10</td>\n",
       "      <td>3.10</td>\n",
       "      <td>3.43</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>0.09</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.08</td>\n",
       "      <td>0.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>0.18</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.15</td>\n",
       "      <td>0.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>0.30</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>1.22</td>\n",
       "      <td>0.09</td>\n",
       "      <td>1.08</td>\n",
       "      <td>1.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>0.13</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.10</td>\n",
       "      <td>0.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>0.24</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.20</td>\n",
       "      <td>0.29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>0.85</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.78</td>\n",
       "      <td>0.94</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>0.41</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.35</td>\n",
       "      <td>0.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>1.01</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.92</td>\n",
       "      <td>1.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>2.09</td>\n",
       "      <td>0.09</td>\n",
       "      <td>1.87</td>\n",
       "      <td>2.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>1.10</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.95</td>\n",
       "      <td>1.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>0.38</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.36</td>\n",
       "      <td>0.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>0.21</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>0.49</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.45</td>\n",
       "      <td>0.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>0.14</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.12</td>\n",
       "      <td>0.16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>1.31</td>\n",
       "      <td>0.06</td>\n",
       "      <td>1.18</td>\n",
       "      <td>1.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>3.88</td>\n",
       "      <td>0.10</td>\n",
       "      <td>3.70</td>\n",
       "      <td>4.04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>0.35</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.31</td>\n",
       "      <td>0.37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>2.56</td>\n",
       "      <td>0.12</td>\n",
       "      <td>2.39</td>\n",
       "      <td>2.79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>53</th>\n",
       "      <td>1.19</td>\n",
       "      <td>0.10</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1.41</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.15</td>\n",
       "      <td>0.27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>55</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>56</th>\n",
       "      <td>0.14</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>57</th>\n",
       "      <td>0.44</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.38</td>\n",
       "      <td>0.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>58</th>\n",
       "      <td>1.42</td>\n",
       "      <td>0.09</td>\n",
       "      <td>1.26</td>\n",
       "      <td>1.56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59</th>\n",
       "      <td>0.55</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.47</td>\n",
       "      <td>0.65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>60</th>\n",
       "      <td>0.46</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.38</td>\n",
       "      <td>0.56</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>61</th>\n",
       "      <td>0.37</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.29</td>\n",
       "      <td>0.44</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62</th>\n",
       "      <td>0.48</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.42</td>\n",
       "      <td>0.58</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>63</th>\n",
       "      <td>0.52</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.46</td>\n",
       "      <td>0.59</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.15</td>\n",
       "      <td>0.23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65</th>\n",
       "      <td>0.32</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>0.59</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.53</td>\n",
       "      <td>0.67</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>68</th>\n",
       "      <td>0.08</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>0.28</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.25</td>\n",
       "      <td>0.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>2.72</td>\n",
       "      <td>0.10</td>\n",
       "      <td>2.60</td>\n",
       "      <td>2.96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>71</th>\n",
       "      <td>1.92</td>\n",
       "      <td>0.09</td>\n",
       "      <td>1.75</td>\n",
       "      <td>2.08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>0.33</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.27</td>\n",
       "      <td>0.39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>73</th>\n",
       "      <td>0.50</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.45</td>\n",
       "      <td>0.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>74</th>\n",
       "      <td>0.09</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75</th>\n",
       "      <td>0.44</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.37</td>\n",
       "      <td>0.59</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>76</th>\n",
       "      <td>0.22</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.20</td>\n",
       "      <td>0.27</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>77</th>\n",
       "      <td>0.11</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.07</td>\n",
       "      <td>0.16</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78</th>\n",
       "      <td>1.39</td>\n",
       "      <td>0.06</td>\n",
       "      <td>1.32</td>\n",
       "      <td>1.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>79</th>\n",
       "      <td>0.36</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>80</th>\n",
       "      <td>0.38</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.32</td>\n",
       "      <td>0.47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>81</th>\n",
       "      <td>0.32</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.29</td>\n",
       "      <td>0.38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82</th>\n",
       "      <td>0.37</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.33</td>\n",
       "      <td>0.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>83</th>\n",
       "      <td>2.00</td>\n",
       "      <td>0.05</td>\n",
       "      <td>1.91</td>\n",
       "      <td>2.12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>84</th>\n",
       "      <td>1.95</td>\n",
       "      <td>0.06</td>\n",
       "      <td>1.82</td>\n",
       "      <td>2.02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>85</th>\n",
       "      <td>0.12</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.09</td>\n",
       "      <td>0.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>86</th>\n",
       "      <td>0.14</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.12</td>\n",
       "      <td>0.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>87</th>\n",
       "      <td>0.30</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.26</td>\n",
       "      <td>0.37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>88</th>\n",
       "      <td>0.19</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.17</td>\n",
       "      <td>0.22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>89</th>\n",
       "      <td>1.48</td>\n",
       "      <td>0.09</td>\n",
       "      <td>1.30</td>\n",
       "      <td>1.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>90</th>\n",
       "      <td>1.22</td>\n",
       "      <td>0.09</td>\n",
       "      <td>1.05</td>\n",
       "      <td>1.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>91</th>\n",
       "      <td>0.11</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.08</td>\n",
       "      <td>0.13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>92</th>\n",
       "      <td>0.41</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.35</td>\n",
       "      <td>0.47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>0.82</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.79</td>\n",
       "      <td>0.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>94</th>\n",
       "      <td>0.54</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.48</td>\n",
       "      <td>0.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>0.34</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.27</td>\n",
       "      <td>0.38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>96</th>\n",
       "      <td>0.29</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>3.37</td>\n",
       "      <td>0.12</td>\n",
       "      <td>3.14</td>\n",
       "      <td>3.53</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>98</th>\n",
       "      <td>0.67</td>\n",
       "      <td>0.08</td>\n",
       "      <td>0.54</td>\n",
       "      <td>0.87</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>99</th>\n",
       "      <td>0.20</td>\n",
       "      <td>0.02</td>\n",
       "      <td>0.16</td>\n",
       "      <td>0.23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>100</th>\n",
       "      <td>1.31</td>\n",
       "      <td>0.07</td>\n",
       "      <td>1.15</td>\n",
       "      <td>1.38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>101</th>\n",
       "      <td>0.28</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.24</td>\n",
       "      <td>0.34</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>102</th>\n",
       "      <td>0.27</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.21</td>\n",
       "      <td>0.32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>103</th>\n",
       "      <td>0.11</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.06</td>\n",
       "      <td>0.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>104</th>\n",
       "      <td>0.59</td>\n",
       "      <td>0.05</td>\n",
       "      <td>0.52</td>\n",
       "      <td>0.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>105</th>\n",
       "      <td>0.41</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.35</td>\n",
       "      <td>0.47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>106</th>\n",
       "      <td>2.22</td>\n",
       "      <td>0.14</td>\n",
       "      <td>2.01</td>\n",
       "      <td>2.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107</th>\n",
       "      <td>0.51</td>\n",
       "      <td>0.04</td>\n",
       "      <td>0.41</td>\n",
       "      <td>0.57</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108</th>\n",
       "      <td>0.43</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.36</td>\n",
       "      <td>0.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>109</th>\n",
       "      <td>0.15</td>\n",
       "      <td>0.03</td>\n",
       "      <td>0.11</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     AVG [%]  ±STD [%]  -CI [%]  +CI [%]\n",
       "0       3.11      0.12     2.87     3.30\n",
       "1       4.66      0.09     4.53     4.79\n",
       "2       0.58      0.05     0.50     0.65\n",
       "3       0.55      0.07     0.47     0.71\n",
       "4       2.54      0.13     2.25     2.70\n",
       "5       0.20      0.02     0.17     0.22\n",
       "6       1.19      0.06     1.10     1.30\n",
       "7       1.57      0.12     1.40     1.81\n",
       "8       4.62      0.11     4.41     4.76\n",
       "9       1.24      0.11     0.98     1.37\n",
       "10      0.78      0.06     0.70     0.92\n",
       "11      0.41      0.04     0.36     0.47\n",
       "12      0.59      0.04     0.52     0.63\n",
       "13      0.19      0.02     0.15     0.21\n",
       "14      2.23      0.11     2.06     2.43\n",
       "15      0.93      0.06     0.84     1.03\n",
       "16      4.62      0.17     4.29     5.01\n",
       "17      0.73      0.04     0.69     0.84\n",
       "18      0.40      0.03     0.34     0.45\n",
       "19      0.68      0.05     0.61     0.78\n",
       "20      0.53      0.03     0.45     0.59\n",
       "21      0.09      0.01     0.07     0.10\n",
       "22      0.55      0.07     0.41     0.66\n",
       "23      2.60      0.12     2.38     2.75\n",
       "24      3.23      0.11     3.07     3.40\n",
       "25      0.08      0.01     0.06     0.10\n",
       "26      1.79      0.06     1.67     1.88\n",
       "27      0.27      0.03     0.21     0.33\n",
       "28      0.62      0.06     0.54     0.77\n",
       "29      0.26      0.02     0.24     0.29\n",
       "30      0.32      0.03     0.26     0.37\n",
       "31      3.32      0.10     3.10     3.43\n",
       "32      0.09      0.01     0.08     0.10\n",
       "33      0.18      0.02     0.15     0.22\n",
       "34      0.30      0.03     0.25     0.36\n",
       "35      1.22      0.09     1.08     1.42\n",
       "36      0.13      0.01     0.10     0.15\n",
       "37      0.24      0.03     0.20     0.29\n",
       "38      0.85      0.04     0.78     0.94\n",
       "39      0.41      0.03     0.35     0.46\n",
       "40      1.01      0.05     0.92     1.14\n",
       "41      2.09      0.09     1.87     2.19\n",
       "42      1.10      0.07     0.95     1.19\n",
       "43      0.38      0.02     0.36     0.42\n",
       "44      0.21      0.03     0.16     0.28\n",
       "45      0.49      0.03     0.45     0.53\n",
       "46      0.14      0.01     0.12     0.16\n",
       "47      1.31      0.06     1.18     1.42\n",
       "48      0.19      0.02     0.16     0.22\n",
       "49      3.88      0.10     3.70     4.04\n",
       "50      0.35      0.02     0.31     0.37\n",
       "51      2.56      0.12     2.39     2.79\n",
       "52      0.20      0.02     0.17     0.23\n",
       "53      1.19      0.10     1.00     1.41\n",
       "54      0.20      0.03     0.15     0.27\n",
       "55      0.19      0.02     0.17     0.23\n",
       "56      0.14      0.02     0.11     0.17\n",
       "57      0.44      0.04     0.38     0.54\n",
       "58      1.42      0.09     1.26     1.56\n",
       "59      0.55      0.05     0.47     0.65\n",
       "60      0.46      0.05     0.38     0.56\n",
       "61      0.37      0.05     0.29     0.44\n",
       "62      0.48      0.05     0.42     0.58\n",
       "63      0.52      0.05     0.46     0.59\n",
       "64      0.20      0.02     0.15     0.23\n",
       "65      0.32      0.05     0.24     0.40\n",
       "66      0.59      0.04     0.53     0.67\n",
       "67      0.19      0.02     0.16     0.23\n",
       "68      0.08      0.01     0.06     0.10\n",
       "69      0.28      0.03     0.25     0.33\n",
       "70      2.72      0.10     2.60     2.96\n",
       "71      1.92      0.09     1.75     2.08\n",
       "72      0.33      0.04     0.27     0.39\n",
       "73      0.50      0.02     0.45     0.53\n",
       "74      0.09      0.01     0.07     0.11\n",
       "75      0.44      0.06     0.37     0.59\n",
       "76      0.22      0.02     0.20     0.27\n",
       "77      0.11      0.02     0.07     0.16\n",
       "78      1.39      0.06     1.32     1.53\n",
       "79      0.36      0.03     0.32     0.40\n",
       "80      0.38      0.04     0.32     0.47\n",
       "81      0.32      0.03     0.29     0.38\n",
       "82      0.37      0.02     0.33     0.42\n",
       "83      2.00      0.05     1.91     2.12\n",
       "84      1.95      0.06     1.82     2.02\n",
       "85      0.12      0.02     0.09     0.15\n",
       "86      0.14      0.02     0.12     0.17\n",
       "87      0.30      0.03     0.26     0.37\n",
       "88      0.19      0.02     0.17     0.22\n",
       "89      1.48      0.09     1.30     1.60\n",
       "90      1.22      0.09     1.05     1.33\n",
       "91      0.11      0.02     0.08     0.13\n",
       "92      0.41      0.03     0.35     0.47\n",
       "93      0.82      0.04     0.79     0.90\n",
       "94      0.54      0.04     0.48     0.60\n",
       "95      0.34      0.04     0.27     0.38\n",
       "96      0.29      0.03     0.24     0.33\n",
       "97      3.37      0.12     3.14     3.53\n",
       "98      0.67      0.08     0.54     0.87\n",
       "99      0.20      0.02     0.16     0.23\n",
       "100     1.31      0.07     1.15     1.38\n",
       "101     0.28      0.03     0.24     0.34\n",
       "102     0.27      0.03     0.21     0.32\n",
       "103     0.11      0.03     0.06     0.17\n",
       "104     0.59      0.05     0.52     0.68\n",
       "105     0.41      0.04     0.35     0.47\n",
       "106     2.22      0.14     2.01     2.46\n",
       "107     0.51      0.04     0.41     0.57\n",
       "108     0.43      0.03     0.36     0.46\n",
       "109     0.15      0.03     0.11     0.20"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1mSettings of this run\u001b[0m\n",
      "Inference mode: peptide\n",
      "EM_convergence_minimum: 0.0001\n",
      "EM_loopcounter_max: 10\n",
      "bootstrap_sampled_fraction: 0.8\n",
      "n_bootstrap_runs: 10\n",
      "CI_percent: 95.0 \n",
      "\n",
      "Runtime: 8.66 seconds ---\n"
     ]
    }
   ],
   "source": [
    "inference_mode = \"peptide\"\n",
    "\n",
    "p_xi_given_zj = np.genfromtxt(\"110 pept 10 prot set/1_p_xi_given_zj_pept_infer_input.csv\", delimiter=',')\n",
    "\n",
    "start_time = time.time() # to start measuring runtime\n",
    "bootstrap_EM()\n",
    "output()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "8c497c21-a39d-40c7-95b1-eebbe3423193",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bootstrap run #0. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #1. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #2. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #3. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #4. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #5. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #6. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #7. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #8. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n",
      "Bootstrap run #9. EM loops: 10\n",
      "WARNING: EM_loopcounter was reached, convergence likely still has not been reached. Consider increasing the maximum number of EM loops.\n"
     ]
    },
    {
     "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>AVG [%]</th>\n",
       "      <th>±STD [%]</th>\n",
       "      <th>-CI [%]</th>\n",
       "      <th>+CI [%]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>28.40</td>\n",
       "      <td>0.35</td>\n",
       "      <td>27.62</td>\n",
       "      <td>28.81</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>24.09</td>\n",
       "      <td>0.45</td>\n",
       "      <td>23.29</td>\n",
       "      <td>24.75</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>15.89</td>\n",
       "      <td>0.28</td>\n",
       "      <td>15.47</td>\n",
       "      <td>16.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7.10</td>\n",
       "      <td>0.37</td>\n",
       "      <td>6.55</td>\n",
       "      <td>7.78</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.93</td>\n",
       "      <td>0.27</td>\n",
       "      <td>5.57</td>\n",
       "      <td>6.38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>7.39</td>\n",
       "      <td>0.45</td>\n",
       "      <td>6.78</td>\n",
       "      <td>8.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3.40</td>\n",
       "      <td>0.42</td>\n",
       "      <td>2.68</td>\n",
       "      <td>4.38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>2.37</td>\n",
       "      <td>0.30</td>\n",
       "      <td>2.05</td>\n",
       "      <td>2.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2.79</td>\n",
       "      <td>0.22</td>\n",
       "      <td>2.49</td>\n",
       "      <td>3.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>2.64</td>\n",
       "      <td>0.28</td>\n",
       "      <td>2.36</td>\n",
       "      <td>3.08</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   AVG [%]  ±STD [%]  -CI [%]  +CI [%]\n",
       "0    28.40      0.35    27.62    28.81\n",
       "1    24.09      0.45    23.29    24.75\n",
       "2    15.89      0.28    15.47    16.30\n",
       "3     7.10      0.37     6.55     7.78\n",
       "4     5.93      0.27     5.57     6.38\n",
       "5     7.39      0.45     6.78     8.10\n",
       "6     3.40      0.42     2.68     4.38\n",
       "7     2.37      0.30     2.05     2.89\n",
       "8     2.79      0.22     2.49     3.19\n",
       "9     2.64      0.28     2.36     3.08"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1mSettings of this run\u001b[0m\n",
      "Inference mode: protein\n",
      "EM_convergence_minimum: 0.0001\n",
      "EM_loopcounter_max: 10\n",
      "bootstrap_sampled_fraction: 0.8\n",
      "n_bootstrap_runs: 10\n",
      "CI_percent: 95.0 \n",
      "\n",
      "Runtime: 0.03 seconds ---\n"
     ]
    }
   ],
   "source": [
    "inference_mode = \"protein\"\n",
    "\n",
    "dyeseq_in_prot_match_count = np.genfromtxt(\"110 pept 10 prot set/2_dyeseq_in_prot_match_count.csv\", delimiter=',') # To calculate this the target proteome needs to be virtually digested and labelled (i.e. turned into dye sequences). Next, the reads' dye sequences are matched against this, and the number of matches per protein is counted. I did not write any code to solve this particular problem, but it seems like it should be fairly easy\n",
    "p_xi_given_zj = np.genfromtxt(\"110 pept 10 prot set/2_p_xi_given_zj_prot_infer_input.csv\", delimiter=',') # The likelihood for any dye seq Xi given a protein Zj is equal to the number of matches of that dye seq in a particular protein, divided by the sum of all dye seq matches of that protein. I did not write any code to solve this problem either\n",
    "\n",
    "start_time = time.time() # to start measuring runtime\n",
    "bootstrap_EM()\n",
    "output()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b5c4ece-e9c8-44c4-987a-c7bd1b15ec6f",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": [],
    "toc-hr-collapsed": true
   },
   "source": [
    "# Legacy stuff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea95e8dc-bea2-4a9e-9927-a24e879c09fe",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Generation of random distributions from 3 Gaussians and plotting them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "b966b932-583d-4442-bd65-6929ad4c8caa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def create_data_from_3_gaussians(mean1, mean2, mean3, std1, std2, std3, n1, n2, n3): # Creation of random data points from multiple Gaussians ki\n",
    "    k1 = np.random.normal(mean1, std1, n1) #creates array with values created through Gaussian\n",
    "    k2 = np.random.normal(mean2, std2, n2)\n",
    "    k3 = np.random.normal(mean3, std3, n3)\n",
    "\n",
    "    return np.concatenate([k1, k2, k3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "298a5777-5f78-470f-8ec6-1c6b31d68e81",
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_xi_given_zj_from_gaussian_datapoints(kall, mean1, mean2, mean3, std1, std2, std3):\n",
    "    # Calculating p_xi_given_zj (this is what Matt is working on with the simulated data)\n",
    "\n",
    "    pdf_probability_k1 = scipy.stats.norm.pdf(kall, loc=mean1, scale=std1)\n",
    "    pdf_probability_k2 = scipy.stats.norm.pdf(kall, loc=mean2, scale=std2)\n",
    "    pdf_probability_k3 = scipy.stats.norm.pdf(kall, loc=mean3, scale=std3)\n",
    "\n",
    "    p_xi_given_zj = np.vstack((pdf_probability_k1,pdf_probability_k2,pdf_probability_k3))\n",
    "    p_xi_given_zj = np.transpose(p_xi_given_zj)\n",
    "\n",
    "    print(\"p_xi_given_zj\")\n",
    "    display(pd.DataFrame(p_xi_given_zj))\n",
    "    print(\"\\n\")\n",
    "    \n",
    "    print(\"likeliest peptide z of each datapoint x according to scipy.stats.norm.pdf (should be more accurate than EM because its dedicated to Gaussians)\")\n",
    "    print(np.argmax(p_xi_given_zj, axis=1)) # for each column of P(Xi|Zj), the most likely Peptide is returned\n",
    "    print(\"\\n\")\n",
    "    \n",
    "    return p_xi_given_zj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "2bd649d2-aa10-4e72-8649-5f290e8bda68",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_3_gaussians_and_calculate_p_xi_given_zj(mean1, mean2, mean3, std1, std2, std3, n1, n2, n3):\n",
    "    kall = create_data_from_3_gaussians(mean1, mean2, mean3, std1, std2, std3, n1, n2, n3)\n",
    "    p_xi_given_zj = p_xi_given_zj_from_gaussian_datapoints(kall, mean1, mean2, mean3, std1, std2, std3)\n",
    "    \n",
    "    return p_xi_given_zj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "fc2bb9da-3d95-4ca4-bf1c-343b61fbfdf4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_histograms_and_pdfs_from_gaussians(mean1, mean2, mean3, std1, std2, std3, n1, n2, n3):\n",
    "    k1 = np.random.normal(mean1, std1, n1) #creates array with values created through Gaussian\n",
    "    k2 = np.random.normal(mean2, std2, n2)\n",
    "    k3 = np.random.normal(mean3, std3, n3)\n",
    "    \n",
    "    # plotting histograms\n",
    "    nbins = 50\n",
    "    plt.hist(k1, label = \"Peptide 0\", bins=nbins, alpha=0.3, density=True, color=\"orange\") # alpha=transparency, density=True normalises to 1 \n",
    "    plt.hist(k2, label = \"Peptide 1\", bins=nbins, alpha=0.3, density=True, color=\"green\")\n",
    "    plt.hist(k3, label = \"Peptide 2\", bins=nbins, alpha=0.3, density=True, color=\"blue\")\n",
    "\n",
    "    # PDF plot\n",
    "    xmin, xmax = plt.xlim() #finds lower and upper bounds of histogram data\n",
    "    x = np.linspace(start=xmin, stop=xmax, num=100) #num is the number of returned data points - the more points, the finer the fit is plotted\n",
    "    p1 = scipy.stats.norm.pdf(x, mean1, std1)\n",
    "    p2 = scipy.stats.norm.pdf(x, mean2, std2)\n",
    "    p3 = scipy.stats.norm.pdf(x, mean3, std3)\n",
    "\n",
    "    plt.plot(x, p1, linewidth=2, color = \"orange\", label = \"Gauss function k1: mean = {:.2f}, STD = {:.2f}\".format(mean1, std1))\n",
    "    plt.plot(x, p2, linewidth=2, color = \"green\", label = \"Gauss function k2: mean = {:.2f}, STD = {:.2f}\".format(mean2, std2))\n",
    "    plt.plot(x, p3, linewidth=2, color = \"blue\", label = \"Gauss function k3: mean = {:.2f}, STD = {:.2f}\".format(mean3, std3))\n",
    "\n",
    "    plt.legend(loc='upper right')\n",
    "    plt.title(\"PDFs of dwarves and humans\")\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddc5c9d6-3030-40ba-bde9-cf052dfce0ec",
   "metadata": {
    "tags": []
   },
   "source": [
    "## EM functions WITHOUT broadcasting (slow as hell!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "2bc33ad4-d4db-4a60-8b68-3a33632c3815",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def NOBROADCAST_update_p_zj_given_xi(p_xi_given_zj, p_zj, p_zj_given_xi):\n",
    "    ## for PEPTIDE inference\n",
    "    if inference_mode == \"peptide\" or \"protein\":\n",
    "        denominator = 0\n",
    "        for i, row in enumerate(p_xi_given_zj): # Calculating/Updating P(Zj|Xi)\n",
    "            #print(\"ROW of P(Xi|Zj):\", i)\n",
    "            for j, cell in enumerate(row):\n",
    "                #print(\"COLUMN of P(Xi|Zj):\", j)\n",
    "                numerator = cell * p_zj[j]\n",
    "                #print(\"numerator:\", numerator)\n",
    "\n",
    "                for l, cell in enumerate(p_zj):\n",
    "                    # print(\"cell\", i, l, p_xi_given_zj[i, l], end=\"\")\n",
    "                    # print(\" * zl\", p_zj[l])\n",
    "                    denominator = denominator + p_xi_given_zj[i, l] * p_zj[l]\n",
    "                \n",
    "                # print(\"p_zj_given_xi[j][i] = numerator/denominator\", numerator, \"/\", denominator)\n",
    "                #display(pd.DataFrame(p_zj_given_xi))\n",
    "                #p_zj_given_xi[j][i] = numerator/denominator\n",
    "                p_zj_given_xi = np.divide(numerator, denominator, out=np.zeros_like(numerator), where=denominator!=0) # if denominator = 0, return 0. This happens sometimes when the likelihood for Zj is equal to 0. This could alternatively be avoided by cropping the rows where the input P(Xi|Zj) scores equal to 0\n",
    "        \n",
    "                \n",
    "                \n",
    "                denominator = 0                                \n",
    "        #display(pd.DataFrame(p_zj_given_xi))    \n",
    "        #print((np.argmax(p_zj_given_xi, axis=1))) # reports index of max value from each row\n",
    "\n",
    "        \n",
    "    ## for PROTEIN inference (denominator is has no loop)\n",
    "    elif inference_mode == \"not protein idk why this doesnt work\":\n",
    "        # print(\"p_xi_given_zj (subarray)\")\n",
    "        # display(pd.DataFrame(p_xi_given_zj))\n",
    "        for i, row in enumerate(p_xi_given_zj): # Calculating/Updating P(Zj|Xi)#\n",
    "            # print(\"ROW of P(Xi|Zj):\", i)\n",
    "            for j, cell in enumerate(row):\n",
    "                # print(\"COLUMN of P(Xi|Zj):\", j)\n",
    "                # print(\"p_zj_given_xi\", cell, \"*\",  p_zj[j], \"/\", p_xi[i],\"=\", cell*p_zj[j]/p_xi[i])\n",
    "                p_zj_given_xi[j][i] = cell * p_zj[j]/p_xi[i]\n",
    "\n",
    "        # display(pd.DataFrame(p_zj_given_xi))\n",
    "       # print((np.argmax(p_zj_given_xi, axis=1))) # reports index of max value from each row\n",
    "        \n",
    "    else:\n",
    "        print(\"ERROR: inference_mode was not set to peptide or protein\")\n",
    "        \n",
    "    return p_zj_given_xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "961d9e06-9f10-4dbf-bbac-6b7ea0ec6f60",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def NOBROADCAST_update_p_zj(p_zj, p_zj_given_xi):  \n",
    "    if inference_mode == \"peptide\":\n",
    "        for j, element in enumerate(p_zj): #updating the expectation value of Zi\n",
    "            p_zj[j] = p_zj_given_xi[j].sum()/p_xi_given_zj.shape[1] # divide by number of peptides z\n",
    "\n",
    "        p_zj_sum = np.sum(p_zj, axis = 0)\n",
    "        p_zj = p_zj/p_zj_sum\n",
    "\n",
    "        # print(\"EM loop\", EM_loopcounter, \")\")\n",
    "        # display(pd.DataFrame(p_zj))\n",
    "    \n",
    "    elif inference_mode == \"protein\":\n",
    "        p_zj.fill(0) # overwrite all values of pzj with 0\n",
    "\n",
    "        print(\"p_zj_given_xi after update\")\n",
    "        display(pd.DataFrame(p_zj_given_xi))\n",
    "        \n",
    "        for j, row in enumerate(p_zj_given_xi):            \n",
    "            for i, element in enumerate(row):\n",
    "                print(\"p_zj_given_xi[j][i] * p_xi[i]\", p_zj_given_xi[j][i], p_xi[i], p_zj_given_xi[j][i] * p_xi[i])\n",
    "                p_zj[j] = p_zj[j] + (p_zj_given_xi[j][i] * p_xi[i])\n",
    "        \n",
    "        p_zj_sum = np.sum(p_zj, axis = 0)\n",
    "        p_zj = p_zj/p_zj_sum\n",
    "        print(\"p_zj (EM loop\", EM_loopcounter, \")\")\n",
    "        display(pd.DataFrame(p_zj))\n",
    "    \n",
    "    else:\n",
    "        print(\"ERROR: inference_mode was not set to peptide or protein\")\n",
    "    \n",
    "    return p_zj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd380528-a6f9-4427-9c4f-fab563fff82c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f84a31f-527a-49da-9cfd-dae2834f9f97",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb3b7723-efb1-4763-9796-904390e1a2f9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}