nbs/02_io.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| default_exp io"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Input / Output\n",
    "\n",
    "AlphaPept can read input from different files. The most important files are those that contain raw MS data. Currently, AlphaPept supports the following raw input formats:\n",
    "\n",
    "* `.raw` files for `Thermo` raw data.\n",
    "* `.d` folders for `Bruker` raw data.\n",
    "* `.wiff` files for `SciEx` raw data.\n",
    "* `.mzML` for generic input.\n",
    "\n",
    "Since 0.5.0, AlphaPept uses [AlphaTims](https://github.com/MannLabs/alphatims) and [AlphaRaw](https://github.com/MannLabs/alpharaw) to import Bruker and Thermo files.\n",
    "\n",
    "After reading these input files, AlphaPept uses several preprocessing functions to prepare it for subsequent analysis. Functions include, for instance:\n",
    "\n",
    "* Centroiding of profile data to centroided data.\n",
    "* Trimming spectra to retain only the `n` most intense peaks.\n",
    "\n",
    "Finally, AlphaPept exports raw data as a `.ms_data.hdf` file. All subsequent results can be added to this file to ensure ease of use/access, full transparency, and reproducibility using a single well-structured file for the whole analysis instead of using randomly scattered files that all relate to a single dataset."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing data\n",
    "\n",
    "Raw data needs to be preprocessed to simplify it for AlphaPept analysis. Currently, this requires:\n",
    "\n",
    "* Profiling of the data\n",
    "* Trimming of spectra to retain only the 'n'-most intense peaks\n",
    "\n",
    "### Centroiding profile data\n",
    "\n",
    "Importing raw data frequently results in profile data. When having profile data, Alphapept first needs to perform centroiding to use this data properly. For this, it needs to search local maxima (\"peaks\") of the intensity as a function of m/z. For this AlphaPept uses the function `get_peaks`. A peak is described by three points, the `start` of the peak, the `center`, and the `end`. The function accepts an intensity array and calculates the delta (gradient) between consecutive data points to determine the start, center, and end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "from numba import njit\n",
    "import numpy as np\n",
    "\n",
    "@njit\n",
    "def get_peaks(int_array: np.ndarray) -> list:\n",
    "    \"\"\"Detects peaks in an array.\n",
    "\n",
    "    Args:\n",
    "        int_array (np.ndarray): An array with intensity values.\n",
    "\n",
    "    Returns:\n",
    "        list: A regular Python list with all peaks.\n",
    "            A peak is a triplet of the form (start, center, end)\n",
    "\n",
    "    \"\"\"\n",
    "    peaklist = []\n",
    "    gradient = np.diff(int_array)\n",
    "    start, center, end = -1, -1, -1\n",
    "\n",
    "    for i in range(len(gradient)):\n",
    "\n",
    "        grad = gradient[i]\n",
    "\n",
    "        if (end == -1) & (center == -1):  # No end and no center yet\n",
    "            if grad <= 0:  # If decreasing, move start point\n",
    "                start = i\n",
    "            else:  # If increasing set as center point\n",
    "                center = i\n",
    "\n",
    "        if (end == -1) & (\n",
    "            center != -1\n",
    "        ):  # If we have a centerpoint and it is still increasing set as new center\n",
    "            if grad >= 0:\n",
    "                center = i\n",
    "            else:  # If it is decreasing set as endpoint\n",
    "                end = i\n",
    "\n",
    "        if end != -1:  # If we have and endpoint and it is going down\n",
    "            if grad < 0:\n",
    "                end = i  # Set new Endpoint\n",
    "            else:  # if it stays the same or goes up set a new peak\n",
    "                peaklist.append((start + 1, center + 1, end + 1))\n",
    "                start, center, end = end, -1, -1  # Reset start, center, end\n",
    "\n",
    "    if end != -1:\n",
    "        peaklist.append((start + 1, center + 1, end + 1))\n",
    "\n",
    "    return peaklist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick test to illustrate the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(2, 6, 9), (9, 13, 16)]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "int_array = np.array([0,0,0,1,2,3,3,2,1,0,0,1,2,3,2,1,0])\n",
    "get_peaks(int_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "def test_get_peaks():\n",
    "    int_array = np.array([0,0,0,1,2,3,3,2,1,0,0,1,2,3,2,1,0])\n",
    "    assert get_peaks(int_array) == [(2, 6, 9), (9, 13, 16)]\n",
    "    \n",
    "test_get_peaks()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine the center of the peak, we distinguish based on the number of raw data points that are contained in the peak:\n",
    "\n",
    "* One data point: m/z value of this point is taken as the center\n",
    "* Two datapoints: Average of the two m/z values, weighted by the intensities\n",
    "* Three or more data points: Gaussian estimation of the center position.\n",
    "\n",
    "For the Gaussian estimation, only the three central points are used to fit a Gaussian Peak shape. The Gaussian is then approximated with the logarithm. \n",
    "\n",
    "The gaussian estimator is defined in `gaussian_estimator` and is used by the wrapper `get_centroid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "from numba import njit\n",
    "\n",
    "@njit\n",
    "def get_centroid(\n",
    "    peak: tuple,\n",
    "    mz_array: np.ndarray,\n",
    "    int_array: np.ndarray\n",
    ") -> tuple:\n",
    "    \"\"\"Wrapper to estimate centroid center positions.\n",
    "\n",
    "    Args:\n",
    "        peak (tuple): A triplet of the form (start, center, end)\n",
    "        mz_array (np.ndarray): An array with mz values.\n",
    "        int_array (np.ndarray): An array with intensity values.\n",
    "\n",
    "    Returns:\n",
    "        tuple: A tuple of the form (center, intensity)\n",
    "    \"\"\"\n",
    "    start, center, end = peak\n",
    "    mz_int = np.sum(int_array[start + 1 : end])\n",
    "\n",
    "    peak_size = end - start - 1\n",
    "\n",
    "    if peak_size == 1:\n",
    "        mz_cent = mz_array[center]\n",
    "    elif peak_size == 2:\n",
    "        mz_cent = (\n",
    "            mz_array[start + 1] * int_array[start + 1]\n",
    "            + mz_array[end - 1] * int_array[end - 1]\n",
    "        ) / (int_array[start + 1] + int_array[end - 1])\n",
    "    else:\n",
    "        mz_cent = gaussian_estimator(peak, mz_array, int_array)\n",
    "\n",
    "    return mz_cent, mz_int\n",
    "\n",
    "@njit\n",
    "def gaussian_estimator(\n",
    "    peak: tuple,\n",
    "    mz_array: np.ndarray,\n",
    "    int_array: np.ndarray\n",
    ") -> float:\n",
    "    \"\"\"Three-point gaussian estimator.\n",
    "\n",
    "    Args:\n",
    "        peak (tuple): A triplet of the form (start, center, end)\n",
    "        mz_array (np.ndarray): An array with mz values.\n",
    "        int_array (np.ndarray): An array with intensity values.\n",
    "\n",
    "    Returns:\n",
    "        float: The gaussian estimate of the center.\n",
    "    \"\"\"\n",
    "    start, center, end = peak\n",
    "\n",
    "    m1, m2, m3 = mz_array[center - 1], mz_array[center], mz_array[center + 1]\n",
    "    i1, i2, i3 = int_array[center - 1], int_array[center], int_array[center + 1]\n",
    "\n",
    "    if i1 == 0:  # Case of sharp flanks\n",
    "        m = (m2 * i2 + m3 * i3) / (i2 + i3)\n",
    "    elif i3 == 0:\n",
    "        m = (m1 * i1 + m2 * i2) / (i1 + i2)\n",
    "    else:\n",
    "        l1, l2, l3 = np.log(i1), np.log(i2), np.log(i3)\n",
    "        m = (\n",
    "            ((l2 - l3) * (m1 ** 2) + (l3 - l1) * (m2 ** 2) + (l1 - l2) * (m3 ** 2))\n",
    "            / ((l2 - l3) * (m1) + (l3 - l1) * (m2) + (l1 - l2) * (m3))\n",
    "            * 1\n",
    "            / 2\n",
    "        )\n",
    "\n",
    "    return m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function performs as expected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.499999999999999\n",
      "(4.499999999999999, 12)\n"
     ]
    }
   ],
   "source": [
    "int_array = np.array([0,0,0,1,2,3,3,2,1,0,0,1,2,3,2,1,0])\n",
    "mz_array = np.array([0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])\n",
    "peak = (2, 6, 9)\n",
    "\n",
    "print(gaussian_estimator(peak, mz_array, int_array))\n",
    "print(get_centroid(peak, mz_array, int_array))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "def test_get_centroid():\n",
    "    \n",
    "    int_array = np.array([0,0,0,1,2,3,3,2,1,0,0,1,2,3,2,1,0])\n",
    "    mz_array = np.array([0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])\n",
    "    \n",
    "    peak = (2, 6, 9)\n",
    "    \n",
    "    mz, intsum = get_centroid(peak, mz_array, int_array)\n",
    "    \n",
    "    assert np.allclose(mz, 4.499999999999999)\n",
    "    assert np.allclose(intsum, 12)\n",
    "    \n",
    "    def gaussian(x, mu, sigma):\n",
    "        return np.exp(-np.power(x - mu, 2.0) / (2 * np.power(sigma, 2.0)))\n",
    "\n",
    "    test_values = np.arange(0, 200)\n",
    "    gauss_trace = np.zeros(len(test_values))\n",
    "\n",
    "    centers = [30, 60, 120]\n",
    "    sigmas = [1, 2, 3]\n",
    "\n",
    "    for center, sigma in zip(centers, sigmas):\n",
    "        gauss_trace += np.round(gaussian(test_values, center, sigma), 3)\n",
    "\n",
    "    peaks = get_peaks(gauss_trace)\n",
    "    centers_found = [_[1] for _ in peaks]\n",
    "\n",
    "    assert set(centers) == set(centers_found)\n",
    "\n",
    "    for i, peak in enumerate(peaks):\n",
    "        center_est, a = get_centroid(peak, test_values, gauss_trace)\n",
    "        sigma_est = (a * 0.3989) / 1\n",
    "\n",
    "        assert np.round(center_est, 2) == centers[i]\n",
    "        assert np.round(sigma_est, 2) == sigmas[i]\n",
    "        \n",
    "test_get_centroid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The detected centroid can also easily be visualized:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAE8CAYAAACRnDftAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAOUlEQVR4nO3deVhTZ9o/8O9JIAn7IrsgqCDKptatYKu4UsVWpjPaetnKWK1OB6daZ8YOnZlW7WuZdl5r+2rHZWyrHcdfbW21rSgWcWlVqriDC4iyCYRFJRCWAMnz+yMkGgGBkHBy4P5c17k0Jyfn3ITk5tznec7zcIwxBkIIERgR3wEQQogxKHkRQgSJkhchRJAoeRFCBImSFyFEkCh5EUIEiZIXIUSQKHkRQgSJkhchRJAoeRHSDo7jsHr16g63W716NTiOM39AxAAlL2ISt27dwtKlSzFo0CDIZDI4Ojpi/Pjx+Pjjj1FfX2+WY5aUlGD16tW4dOmSWfZPLJsV3wEQ4UtOTsacOXMglUqxYMEChIWFobGxESdPnsSf//xnXL16Fdu2bTP5cUtKSrBmzRoEBARgxIgRJt9/fX09rKzoK2Kp6DdDuiUvLw8vvvgi/P39cfToUXh7e+ufS0hIQG5uLpKTk3mM8IG6ujrY2tp2enuZTGbGaEh3UdlIuuWDDz6AUqnEp59+apC4dAIDA7F8+XL94127dmHUqFGwsbGBq6srXnzxRRQVFRm8Jjo6GmFhYbh27RomTZoEW1tb9O/fHx988IF+m+PHj2PMmDEAgIULF4LjOHAchx07dhjs4/z585gwYQJsbW3x1ltvAQDKy8uxaNEieHp6QiaTYfjw4di5c2er2Nu65nXy5EmMGTMGMpkMgwcPxtatW41630j30ZkX6ZYffvgBgwYNQlRUVIfbrlu3Dn//+98xd+5cLF68GBUVFdi4cSMmTJiAixcvwtnZWb/t/fv38cwzz+D555/H3LlzsXfvXrz55psIDw/HjBkzMGzYMKxduxZvv/02lixZgqeffhoADOK4e/cuZsyYgRdffBEvvfQSPD09UV9fj+joaOTm5mLZsmUYOHAgvv76a/z2t79FVVWVQaJ9VGZmJqZPnw53d3esXr0azc3NeOedd+Dp6Wn8G0iMxwgxkkKhYADY7NmzO9w2Pz+ficVitm7dOoP1mZmZzMrKymD9xIkTGQD2xRdf6NepVCrm5eXFfv3rX+vXZWRkMADs888/b3U83T62bNlisP6jjz5iANiuXbv06xobG1lkZCSzt7dn1dXV+vUA2DvvvKN/HBcXx2QyGSsoKNCvu3btGhOLxYy+Sj2PykZitOrqagCAg4NDh9t+++230Gg0mDt3LiorK/WLl5cXgoKCcOzYMYPt7e3t8dJLL+kfSyQSjB07Frdv3+50fFKpFAsXLjRYd/DgQXh5eWHevHn6ddbW1nj99dehVCpx4sSJNvelVqtx+PBhxMXFYcCAAfr1w4YNQ0xMTKdjIqZDZSMxmqOjIwCgpqamw21v3rwJxhiCgoLafN7a2trgsa+vb6u+Uy4uLrhy5Uqn4+vfvz8kEonBuoKCAgQFBUEkMvy7PWzYMP3zbamoqEB9fX2b8QcHB+PgwYOdjouYBiUvYjRHR0f4+PggKyurw201Gg04jsOhQ4cgFotbPW9vb2/wuK1tAIB1YdRyGxubTm9LhIeSF+mWWbNmYdu2bUhPT0dkZGS72w0ePBiMMQwcOBBDhgwxybGN6dXu7++PK1euQKPRGJx93bhxQ/98W9zd3WFjY4ObN2+2ei47O7vLcZDuo2tepFtWrVoFOzs7LF68GGVlZa2ev3XrFj7++GM8//zzEIvFWLNmTauzJ8YY7t692+Vj29nZAQCqqqo6/ZqZM2dCLpdjz549+nXNzc3YuHEj7O3tMXHixDZfJxaLERMTg/3796OwsFC//vr16zh8+HCXYyfdR2depFsGDx6M3bt344UXXsCwYcMMetifPn1a3w1h+fLl+J//+R8kJiYiPz8fcXFxcHBwQF5eHvbt24clS5bgT3/6U5eP7ezsjC1btsDBwQF2dnYYN24cBg4c2O5rlixZgq1bt+K3v/0tzp8/j4CAAOzduxenTp3CRx999NjGhzVr1iAlJQVPP/00fv/73+uTXmhoaJeuxRET4bWtk/QaOTk57NVXX2UBAQFMIpEwBwcHNn78eLZx40bW0NCg3+6bb75hTz31FLOzs2N2dnZs6NChLCEhgWVnZ+u3mThxIgsNDW11jPj4eObv72+w7rvvvmMhISHMysrKoNtEe/tgjLGysjK2cOFC5ubmxiQSCQsPD2+zuwUe6SrBGGMnTpxgo0aNYhKJhA0aNIht2bKFvfPOO9RVggccYzRvIyFEeOiaFyFEkCh5EUIEiZIXIUSQKHkRQgSJkhchRJAoeRFCBKnPdVLVaDQoKSmBg4MDTZpAiAVijKGmpgY+Pj6tbqB/WJ9LXiUlJfDz8+M7DEJIB4qKiuDr69vu830ueelu/ygqKtIP6dLjamsBHx/t/0tKgJZ79EjX1TbWwme99r0s+WMJ7CT0XhrFgj6T1dXV8PPz63CcuD6XvHSloqOjI3/J6+HhXhwdKXl1g7hRDLTMk+Ho6EjJy1gW+Jns6LIOXbAnhAgSJS9CiCBR8iKECBKvyWvz5s2IiIjQX3+KjIzEoUOHHvuar7/+GkOHDoVMJkN4eDiNHU5IH8Vr8vL19cU//vEPnD9/HufOncPkyZMxe/ZsXL16tc3tT58+jXnz5mHRokW4ePEi4uLiEBcX16kx1Ilw3LkDHDum/VdI+yY9jOfxxFpxcXFh27dvb/O5uXPnstjYWIN148aNY0uXLu30/nVzDSoUim7F2S1KJWOAdlEq+YvDAm3fzphIpH1rRCLt48dRqpQMq8GwGkypevx72dV99ykW9Jns7HfUYrpKqNVqfP3116itrW13Iof09HSsXLnSYJ1uXPH2qFQqqFQq/WPdXIPE8ty5AyxZwqDRaJvINRpg8asM/yk6A5lzY5uvaWb1+v/P/uQkrLi2ZwxqqJLgxNpxAHuw76VLgZgY4DH9IIkF4z15ZWZmIjIyEg0NDbC3t8e+ffsQEhLS5rZyubzV1Oqenp6Qy+Xt7j8pKQlr1qwxaczEPG7ehD5x6TEON3IYZAPanhtSgwagJV/llCkhQnOb2zUU9NMnLh21GsjNpeQlVLwnr+DgYFy6dAkKhQJ79+5FfHw8Tpw40W4C66rExESDszVd711ieYKCAHDMIMmIRAz/ei0YHl6aNl9T31yL6V9p//9p/GjYWLXdubJcLsKcr5hBchSLgcBAk4VPehjvyUsikSCw5RM0atQoZGRk4OOPP8bWrVtbbevl5dVqeq2ysjJ4eXm1u3+pVAqpVGraoIlZSJwa0O+ZHNxNCQOYCGIxsHUrh+efdm33NbWND8rEqMFu7fewDwK2bQOWLmVQqzmA0+C99U3w9aXPhlBZXD8vjUZjcI3qYZGRkUhLSzNYl5qa+tjJTolwpGTJYR9RhGlrz+PYMSA/H1i0yHT7X7QIyM/nMOVP19D/d8dgP7zIdDsnPY7X5JWYmIiffvoJ+fn5yMzMRGJiIo4fP4758+cDABYsWIDExET99suXL0dKSgrWr1+PGzduYPXq1Th37hyWLVvG149ATCj5SikAYM6EfoiONs+1KF9f4JXf2MPKsQEHM0tNfwDSY3hNXuXl5ViwYAGCg4MxZcoUZGRk4PDhw5g2bRoAoLCwEKWlDz5gUVFR2L17N7Zt24bhw4dj79692L9/P8LCwvj6EYiJlFc34Gz+PQDAjPD2LwOYQkyoJ8QiDlnF1civrDXrsYj58HrN69NPP33s88ePH2+1bs6cOZgzZ46ZIiJ8OZQlB2PAyAHO8HWxNeux+tlLETW4H36+WYnkzFIkTKKr9kJkcde8SN+kKxljw7175Hi64+iOS4SHkhfhnVzRgIwCbck4s4eSV0yoF8QiDtdKq3G7QtkjxySmRcmL8O5QVikYA0b5u8DHue0e8qbmYifB+EA3AKAL9wJFyYvwrqdLRp1ZLcc7QKWjIFHyIrwqVdTjXMF9AOZvZXzU9FBPWIk43JDXILecSkehoeRFeHUwU3tf6mh/F3g79UzJqONsS6WjkFHyIrzSJY3YiJ4tGXV0x6XkJTyUvAhvSqrqcb7gPjgOmBHGT/KKCfGCtVhXOrY9cgWxTJS8CG90Zztj/F3h5STjJQYnW2s81VI6Jl9pf2glYnkoeRHeJPNcMurERmgnW03OLOE1DtI1lLwIL+7cr8PFwqqWkrFnWxkfNS3EE9ZiDjllSuSUUekoFJS8CC8OtbQyjg1whYcjPyWjjpONNSYEuQOg24WEhJIX4cWBlpJxFs8lo46udE3OLAVjjOdoSGdQ8iI9ruheHS4XaUvGGJ5LRp2pIZ6QiEXILVcip4w6rAoBJS/S4w5lac+6xg10hYcDvyWjjqPMGhOGtLQ6Up8vQaDkRXqc/l7GllY+S6EvHa+UUOkoAJS8SI8quleHy3cUEHHAM6GWUTLqTB3mCYmVCLcqapFNrY4Wj5IX6VG6kuzJQf3g7mBZM/c4yKwxcQi1OgoFJS/Sox6UjJbRyvioWREPRlil0tGyUfIiPabgbi0yiy2zZNSZ0lI63q6sxfVSKh0tGSUv0mN0JWPUYDf0s7esklHHXmqFScEtpSPdLmTRKHmRHmPpJaOO/l5HKh0tGiUv0iPyK2txtaQaYhGHGAstGXWmDPWA1EqE/Lt1uFpSzXc4pB2UvEiPeFAy9oOrnYTnaB7PTmqFScEeAGiQQkvGa/JKSkrCmDFj4ODgAA8PD8TFxSE7O/uxr9mxYwc4jjNYZDLL6KVN2sfXJBvGonsdLR+vyevEiRNISEjAL7/8gtTUVDQ1NWH69OmorX38FOyOjo4oLS3VLwUFBT0UMTHG7QolrpUKo2TUmTzUAzJrEQqodLRYVnwePCUlxeDxjh074OHhgfPnz2PChAntvo7jOHh5CeNLQB6UXuMD3eBi4SWjjp3UCpOHeuBgphwHrpQirL8T3yGRR1jUNS+FQgEAcHV1fex2SqUS/v7+8PPzw+zZs3H16tV2t1WpVKiurjZYSM/SzYs4SyAlo05s+IMRVql0tDwWk7w0Gg1WrFiB8ePHIywsrN3tgoOD8dlnn+G7777Drl27oNFoEBUVhTt37rS5fVJSEpycnPSLn5+fuX4E0obcciVuyGtgJeIwPdST73C6ZNJQd9hYi1F0rx6ZxQq+wyGPsJjklZCQgKysLHz55ZeP3S4yMhILFizAiBEjMHHiRHz77bdwd3fH1q1b29w+MTERCoVCvxQVFZkjfNIOXcn4VJAbnG2FUTLq2EqsMHmYttWR7nW0PBaRvJYtW4YDBw7g2LFj8PX17dJrra2tMXLkSOTm5rb5vFQqhaOjo8FCeo7uSz9TYCWjjq519AB1WLU4vCYvxhiWLVuGffv24ejRoxg4cGCX96FWq5GZmQlvb2F+OXqz3PIaZJfVwFrMISZEmA0sk4I9YGMtRnFVPS7fodLRkvCavBISErBr1y7s3r0bDg4OkMvlkMvlqK+v12+zYMECJCYm6h+vXbsWP/74I27fvo0LFy7gpZdeQkFBARYvXszHj0AeQzcP4lOBbnCyteY5GuPYSMSYMow6rFoiXpPX5s2boVAoEB0dDW9vb/2yZ88e/TaFhYUoLX3wobl//z5effVVDBs2DDNnzkR1dTVOnz6NkJAQPn4E8hi6G5stbcTUrqJhciwTr/28OvNBOH78uMHjDRs2YMOGDWaKiJhKTlkNcsqUsBZzmBYirFbGR0UHe8BWoi0dLxVVYeQAF75DIrCQC/ak99FdqJ8Q5A4nG2GWjDoyazGmDtMmYGp1tByUvIjJMcb0N2Jb+vA3naX7OQ5mlkKjodLRElDyIiaXU6ZEbrkSErEIUwVeMupMHOIOO4kYJYoGXCyq4jscAkpexAySr2gv1E8Y4g5HmbBLRh2ZtVh/7Y5KR8tAyYuYFGMMB/QlozD7drVH19GWSkfLQMmLmNQNeQ1uV9RCYiXSX+TuLSYMcYe91Ary6gZcKLzPdzh9HiUvYlK6jpwTh7jDoZeUjDoGpSN1WOUdJS9iMowx/fWgWb2klfFRsVQ6WgxKXsRkrpfW4HaltmSc0stKRp2nh7jBQWqFsmoVzlPpyCtKXsRkdLcDTQrWXhvqjaRWYkwLpVZHS0DJi5jEwyWj0O9l7Mishzqsqql05A0lL2ISV0uqkX+3DlIrEaYM9eA7HLN6KtAdDjIrlNeocC7/Ht/h9FmUvIhJ6FrfJg/1gF0vLRl1JFYi/SxI1OrIH0pepNseLhmFOmJqVz1odZRT6cgTSl6k266WVKPwXh1k1iJM7uUlo874QDc4yqxQqVQhg0pHXlDyIt2mm9qsL5SMOgalI7U68oKSF+kW7fA3LSOmhvfuVsZH6YbJOZRFrY58oORFuiWzWIGie/WwsRZj0lB3vsPpUeMD3eBkY41KZSPO5N3lO5w+h5IX6RZdyTR5mAdsJX2jZNSxFovwDJWOvKHkRYzGGNNf75rVR1oZH6UrHVOy5GhWa3iOpm+h5EWMdvmOAsVV9bCViBEd3DdaGR8VObgfXGytcbe2EWfyqNWxJ1HyIkbTjZg6ZZgnbCRinqPhh7VYhGfCtKXjASodexQlL2IUxhgOZmonlY0N710jpnaVrmPu4atUOvYkSl7EKJeKqvp8yagTOUhbOt6rbcQvt6l07Cm8Jq+kpCSMGTMGDg4O8PDwQFxcHLKzszt83ddff42hQ4dCJpMhPDwcBw8e7IFoycN0rWtTh3lCZt03S0YdK7EIz4S1zKrd0ueNmB+vyevEiRNISEjAL7/8gtTUVDQ1NWH69Omora1t9zWnT5/GvHnzsGjRIly8eBFxcXGIi4tDVlZWD0betxUWMvy//Q1orpb1mnkZu0s3TM73p+4j9YgGd+7wHFBfwCxIeXk5A8BOnDjR7jZz585lsbGxBuvGjRvHli5d2qljKBQKBoApFIpuxdotSiVjgHZRKvmLwwjbtzMmEmm04XMatnmrmtd4lColw2owrAZTqvh7L5ua1WzA7CwGTvveiETa90owLOgz2dnvqEVd81IoFAAAV1fXdrdJT0/H1KlTDdbFxMQgPT29ze1VKhWqq6sNFmKcO3eAJUsAjYbTrmAclv1eRGcZAOSlIhR9HwIw7Xuj0QBLl4LeGzOymOSl0WiwYsUKjB8/HmFhYe1uJ5fL4elpOD66p6cn5HJ5m9snJSXByclJv/j5+Zk07r7k5k3tl/JhajWQm8tPPJbk5k2AtSQuHXpvzMtikldCQgKysrLw5ZdfmnS/iYmJUCgU+qWoqMik++9LgoIAkcjwBmSxGAgM5CkgC0LvTc+ziOS1bNkyHDhwAMeOHYOvr+9jt/Xy8kJZWZnBurKyMnh5td3XSCqVwtHR0WAhxvH1BWJ/XwJw2tMvsRjYulW7vq/z9QW2bePAtSQwTsTovTEzXpMXYwzLli3Dvn37cPToUQwcOLDD10RGRiItLc1gXWpqKiIjI80VJmmh0TCUeV9H/98dw/9+fg/5+cCiRXxHZTkWLQL2naiC57x0BC8/gZfjqcOqOfGavBISErBr1y7s3r0bDg4OkMvlkMvlqK+v12+zYMECJCYm6h8vX74cKSkpWL9+PW7cuIHVq1fj3LlzWLZsGR8/Qp9yvvA+yqpVcHFvxrKXnOisog2zopzhG1KLekktTt2q5DucXo3X5LV582YoFApER0fD29tbv+zZs0e/TWFhIUpLH9wzFhUVhd27d2Pbtm0YPnw49u7di/379z/2Ij8xDV3H1GmhnpBa9e2Oqe0RizjMDKdhcnoCrwMwMdbx6JPHjx9vtW7OnDmYM2eOGSIi7VFrGA62zJQzizqmPlZsuDe+SC/A4atyvPercEisLOLScq9D7yrplHP591Beo4KDzApPBfatEVO7anSAKzwcpKhpaMbJ3Aq+w+m1KHmRTtHNTxgT6kVnEh3Qlo7as1MaJsd8jPoUxsfH46effjJ1LMRCaUvGluFvqGTsFN37lHq1DKpmNc/R9E5GJS+FQoGpU6ciKCgI7733HoqLi00dF7EgZ/PuoVKpgpONNcYPduM7HEEYNcAFno5S1Kia8XMOtTqag1HJa//+/SguLsZrr72GPXv2ICAgADNmzMDevXvR1NRk6hgJz3TDvEwP8aSSsZNED5WOupKbmJbRn0R3d3esXLkSly9fxpkzZxAYGIiXX34ZPj4+eOONN3Dz5k1Txkl4otYwpGRRyWgMXats6rUyNDRR6Whq3f4zWlpaitTUVKSmpkIsFmPmzJnIzMxESEgINmzYYIoYCY/O5N1FpbJRWzIGUsnYFSP9XODlKINS1Yyfb1LpaGpGJa+mpiZ88803mDVrFvz9/fH1119jxYoVKCkpwc6dO3HkyBF89dVXWLt2ranjJT1M19HymVAvWIupZOwKg9LxCo2wampGdVL19vaGRqPBvHnzcPbsWYwYMaLVNpMmTYKzs3M3wyN8alZrqGTsptgIb3x2Kk9fOvb1IbNNyajktWHDBsyZMwcymazdbZydnZGXl2d0YIR/Z/Lu4W5tI1xsrRE5uB/f4QjSSD9n+DjJUKJowImcCsSE9u2ZlkzJqDrg2LFjbbYq1tbW4pVXXul2UMQy6DpYPhNGJaOxDEtHanU0JaM+kTt37jQY+UGnvr4eX3zxRbeDIvzTlozaL1tsuA/P0QibruQ+cp1aHU2pS2VjdXU1GGNgjKGmpsagbFSr1Th48CA8PPr2HH69Rfrtu7hf1wRXOwmeHNT+nAKkYyP8nNHf2QbFVfU4nl2unyaNdE+XkpezszM4jgPHcRgyZEir5zmOw5o1a0wWHOGPrsSJCfWCFZWM3cJxHGIjvLHtp9s4cKWUkpeJdCl5HTt2DIwxTJ48Gd98843BLD8SiQT+/v7w8aESQ+ia1BocvqptZaThb0wjNlybvI7eKEd9oxo2Emp17K4uJa+JEycCAPLy8jBgwABwHNfBK4gQpd/Sloz97CQYN5BKRlOI8HWCr4sN7tzXlo4zwumPQnd1OnlduXIFYWFhEIlEUCgUyMzMbHfbiIgIkwRH+JH8UCsjlYymwXEcYsO9sfWn2ziQWUrJywQ6nbxGjBgBuVwODw8PjBgxAhzHtTkSKsdxUKupRUWomtQapFyljqnmEBuhTV5Hr5ejrrEZthJeBzIWvE6/e3l5eXB3d9f/n/ROp3Iroahvgpu9BOMGUsdUUwrv7wQ/VxsU3avHsRsV9MehmzqdvPz9/dv8P+lddCXjjDBviEV0TdOUtKWjD7acuIXkzBJKXt1kdCfV5ORk/eNVq1bB2dkZUVFRKCgoMFlwpGc1Nj9oZaQvlnnoWm+P3ihHraqZ52iEzajk9d5778HGxgYAkJ6ejk2bNuGDDz6Am5sb3njjDZMGSHrOqdxKVDc0w91BijEB1MpoDqE+jvDvZ4uGJg2O3ijnOxxBMyp5FRUVITAwEIB2VNXf/OY3WLJkCZKSkvDzzz+bNEDScw7oS0YvKhnNRNfqCNC9jt1lVPKyt7fH3bt3AQA//vgjpk2bBgCQyWRt3vNILF9jswY/XmspGakZ36x0JfmxbCodu8Oo5DVt2jQsXrwYixcvRk5ODmbOnAkAuHr1KgICAjq9n59++gnPPvssfHx8wHEc9u/f/9jtjx8/rr896eFFLpcb82OQh5zMrUBNQzM8HKQYTSWjWYV4O2Kgmx1UzRqkUeloNKOS1yeffILIyEhUVFTgm2++Qb9+2ib18+fPY968eZ3eT21tLYYPH45PPvmkS8fPzs5GaWmpfqGbwbtPVzLODKdWRnMzLB1phFVjGdVLztnZGZs2bWq1vqs3Zc+YMQMzZszo8vE9PDxolFYTUjWrkXq1DAC1MvaUmeHe2HQsF8eyK6BUNcNeSh1Wu8rod6yqqgpnz55FeXk5NBqNfj3HcXj55ZdNElx7RowYAZVKhbCwMKxevRrjx49vd1uVSgWVSqV/XF1dbdbYhOjnnErUqJrh6SjFqAEufIfTJwzzdsAgNzvcrqxF2vUyzB7Rn++QBMeo5PXDDz9g/vz5UCqVcHR0NLhB25zJy9vbG1u2bMHo0aOhUqmwfft2REdH48yZM3jiiSfafE1SUhIN09MB3byCM8O9IaKSsUfohsnZeDQXB66UUvIyglHXvP74xz/ilVdegVKpRFVVFe7fv69f7t27Z+oY9YKDg7F06VKMGjUKUVFR+OyzzxAVFfXYKdYSExOhUCj0S1FRkdniE6KGJjVSr2lLRhr+pmfpSvQT2RWoaaDJmrvKqORVXFyM119/Hba2tqaOp8vGjh2L3Nzcdp+XSqVwdHQ0WMgDP+Vor7l4O8kw0o9Kxp4U7OmAwe52aFRrcOR6Gd/hCI5RySsmJgbnzp0zdSxGuXTpEry96YzBWLqScUYYlYw9TVs6agfvpA6rXWfUNa/Y2Fj8+c9/xrVr1xAeHg5ra2uD55977rlO7UepVBqcNeXl5eHSpUtwdXXFgAEDkJiYiOLiYv2kHh999BEGDhyI0NBQNDQ0YPv27Th69Ch+/PFHY36MPq+hSY0j16iVkU+zIrzxf2k38VNOJaobmuAos+74RQSAkcnr1VdfBYA2Z8Tuynhe586dw6RJk/SPV65cCQCIj4/Hjh07UFpaisLCQv3zjY2N+OMf/4ji4mLY2toiIiICR44cMdgH6bwTORWobVTDx0mGkX7OfIfTJw3xdECQhz1ulitx5FoZnn/Cl++QBMOo5PVw14juiI6ObnNAQ50dO3YYPF61ahVWrVplkmOTB6UKtTLyKzbCGx8duYnkK6WUvLqg22P8NjQ0mCIO0sMamtT6i8RUMvJL19v+p5sVUNRTq2NnGZW81Go13n33XfTv3x/29va4ffs2AODvf/87Pv30U5MGSMzjeHY56hrV6O9sgxFUMvIqyNMBQzzt0aRm+m4rpGNGJa9169Zhx44d+OCDDyCRSPTrw8LCsH37dpMFR8xHdy9jbIQ3zQJlAXSzktO9jp1nVPL64osvsG3bNsyfPx9i8YP554YPH44bN26YLDhiHvWNaqRd145mQMPfWIbYCC8AwM83K6Goo9KxM4zupKobjPBhGo0GTU30xlu6Y9nlqG9Sw9fFBhG+TnyHQwAEejhgqJcDmjUMh6/REE+dYVTyCgkJaXPE1L1792LkyJHdDoqYl66VMTacSkZLQiOsdo1RXSXefvttxMfHo7i4GBqNBt9++y2ys7PxxRdf4MCBA6aOkZhQXWOzfux0amW0LDMjvLE+NQencitRVdcIZ1tJxy/qw4w685o9ezZ++OEHHDlyBHZ2dnj77bdx/fp1/PDDD/ohoYllOnajAvVNavi52iC8P5WMlmSwuz2GeTuiWcPw41VqdeyI0eN5Pf3000hNTTVlLKQHJGdqW7Niw32oZLRAsyK8cb20GgcySzF3jB/f4Vg0o868Bg0apJ+A42FVVVUYNGhQt4Mi5lGrelAy0vA3lmlmy3WvU7mVuF/byHM0ls2o5JWfn9/m/YsqlQrFxcXdDoqYx9Eb5Who0sC/ny1CfWhoIEs00M0OoT6OUGuYfgJg0rYulY3ff/+9/v+HDx+Gk9ODayZqtRppaWldmj2I9CxqZRSGmeHeuFpSjeTMUrw4dgDf4VisLiWvuLg4ANqRI+Lj4w2es7a2RkBAANavX2+y4IjpKFXNOJZNrYxCEBvujX8ezsbpW3dxV6lCP3sp3yFZpC6VjRqNBhqNBgMGDNBPvKFbVCoVsrOzMWvWLHPFSroh7XoZVM0aBPSzRYg3lYyWLMDNDmH9daUjtTq2x6hrXnl5eXBzczN1LMSMDmbSvYxCorvXUfd7I60Z3VUiLS0NaWlpraY+A4DPPvus24ER09GWjBUAHnwpiGWLDffG+yk3cPpWJZWO7TDqzGvNmjWYPn060tLSUFlZaTB70P37900dI+mmtOtlaGzWYJCbHYZ5O/AdDumEAf1sEeHrBA0DUqjVsU1GnXlt2bIFO3bsMPvkssQ0aPgbYYoN98aVOwokXynF/HH+fIdjcYw682psbERUVJSpYyFmUNPQhBO6kpFaGQVF12H1l9t3UVGj6mDrvseo5LV48WLs3r3b1LEQMzhyvQyNag0Gu9sh2JNKRiHxc7XFcD9nKh3bYVTZ2NDQgG3btuHIkSOIiIhoNfXZhx9+aJLgSPfpO6ZG0L2MQjQr3BuXi6qQfKUELz9JpePDjEpeV65cwYgRIwAAWVlZpoyHmJCivgk/5VQCoHsZhWpGuBfWHbyOM3n3UF7TAA8HGd8hWQyjktexY8dMHQcxgyPXtCVjoIc9hlDJKEi+LrYY4eeMS0VVSMmSY0FkAN8hWYwuJa/nn3++w204jsM333xjdEDEdPQdU2mcekGbFeGNS0VVSL5SSsnrIV26YO/k5NTh4ujY+VtPfvrpJzz77LPw8dFej9m/f3+Hrzl+/DieeOIJSKVSBAYGtpqYlmgp6pvw001qZewNZrT88Tmbfw/l1TRPqk6Xzrw+//xzkx68trYWw4cPxyuvvNKps7q8vDzExsbid7/7Hf773/8iLS0Nixcvhre3N2JiYkwam9B9eawSNbddMWQIo5JR4Po72+CJAc44m9WAD3dWYdmvvOBLE2sDzEIAYPv27XvsNqtWrWKhoaEG61544QUWExPT6eMoFAoGgCkUCmPCNA2lkjFAuyiVJt/99u2McZyGAYxxIg3bvt3kh7AYSpWSYTUYVoMpVaZ/Ly1F/JvlDC2/U5GImf53aubPZFd09jtqVD8vvqSnp2Pq1KkG62JiYpCent7ua1QqFaqrqw2W3uzOHWDJEgbGtN0imIbD0qXa9USY7twB/vNPN6Dld6rRgH6nMLKTKl/kcjk8PT0N1nl6eqK6uhr19fVtviYpKcngmpyfX+8eF/zmTUCjMezPpVYDubk8BUS6jX6nbRNU8jJGYmIiFAqFfikqKuI7JLMKCgI4jhmsE4uBNuYIJgIRFASIHvmm0u9UYMnLy8sLZWWGg7OVlZXB0dERNjY2bb5GKpXC0dHRYOnN7Fwa4TYjE+C0wxSJxcDWraALvALm6wts2waIxS1/lDgN3t/Q2Od/p0aP58WHyMhIHDx40GBdamoqIiMjeYrI8vx4TQ7b8CJMHFOH1dFPIjCQEldvsGgREBPDYf6HWbjVUAbnkYMADOQ7LF7xeualVCpx6dIlXLp0CYC2K8SlS5dQWFgIQFvyLViwQL/97373O9y+fRurVq3CjRs38K9//QtfffUV3njjDT7Ct0jJmdobeH8zoR+ioylx9Sa+vsDCX9vByrGBRlgFz8nr3LlzGDlyJEaOHAkAWLlyJUaOHIm3334bAFBaWqpPZAAwcOBAJCcnIzU1FcOHD8f69euxfft26uPV4n5tI07lau9lnEm96nulGeHe4DjgfMF9lFS13UjVV/BaNkZHR4Mx1u7zbfWej46OxsWLF80YlXAdviqHWsMQ4u2IQe72fIdDzMDTUYYxAa44m3cPBzNLsfjpvjvJs6Au2JPHS35okg3Se+lGCEnu46UjJa9e4q5ShdO37gKgG7F7u2fCvMBxwMXCKty5X8d3OLyh5NVLHL5aBrWGIay/IwLc7PgOh5iRh4MM4wa6AgAOZfbdEVYpefUSyZklAGhqs74iNkL7ez7Qh0tHSl69QKVShXQqGfuUZ0K9IOKAy0VVKLrXN0tHSl69QEqWHBoGhPd3woB+tnyHQ3qAu4MU4wb2A9B3Z9Wm5NULHKRWxj5J9/um5EUEqaJGhV9uU8nYFz0T1lI63lH0ydKRkpfApVzVlozDfZ3g50olY1/iZi9F5GBt6dgX+3xR8hK45CstrYxUMvZJutZl3fycfQklLwErr2nAmbx7AOhexr4qJtQTYhGHzGIFCu7W8h1Oj6LkJWApWXIwBozwc4avC5WMfVE/eymi+mjpSMlLwA60lAo0G3bfpmuo6WulIyUvgSqrbkBGvrZknEElY582PdQLYhGHqyXVyKvsO6UjJS+BOpRZCsaAkQOc0d+57SGwSd/gaifRl459qc8XJS+BOthyQy717SLAQ8Pk9KHSkZKXAMkVDcgooFZG8sD0EC9YiThcK63G7Qol3+H0CEpeAnQoS1syjvJ3gQ+VjASAi50E4wPdAPSd0pGSlwDpSgMqGcnDdB2VD/SR0pGSl8CUKupxruA+ACoZiaGYEC9YiznckNcgt7z3l46UvARGd6F+TIALvJxkPEdDLImTrTWe6kOlIyUvgdHfy0hnXaQNuhFW+0KrIyUvASmpqseFwipwHHVMJW2bFuIJazGH7LIa5JbX8B2OWVHyEhBdKTDG3xWejlQyktacbKzxdJA7ACD5Su+enMMiktcnn3yCgIAAyGQyjBs3DmfPnm132x07doDjOINFJusbX2Sal5F0hv5ex5ZJWXor3pPXnj17sHLlSrzzzju4cOEChg8fjpiYGJSXl7f7GkdHR5SWluqXgoKCHoyYH3fu1+GirmQM8+I7HGLBpoZ4QiIWIadMiZyy3ls68p68PvzwQ7z66qtYuHAhQkJCsGXLFtja2uKzzz5r9zUcx8HLy0u/eHp69mDE/NDNzzc2wBUeVDKSx3CyscaEIdpWx9584Z7X5NXY2Ijz589j6tSp+nUikQhTp05Fenp6u69TKpXw9/eHn58fZs+ejatXr7a7rUqlQnV1tcEiRLr5+Wj4G9IZuksLyZmlYIzxHI158Jq8KisroVarW505eXp6Qi5v+2JjcHAwPvvsM3z33XfYtWsXNBoNoqKicOfOnTa3T0pKgpOTk37x8/Mz+c9hbkX36nC5qAoiDoihkpF0wtRhnpBYiZBbrkROWe/ssMp72dhVkZGRWLBgAUaMGIGJEyfi22+/hbu7O7Zu3drm9omJiVAoFPqlqKiohyPuPl0r47iB/eDhQCUj6ZiDzBoTh+haHXvnhXtek5ebmxvEYjHKysoM1peVlcHLq3NnGNbW1hg5ciRyc3PbfF4qlcLR0dFgERpqZSTG0F1iONBLS0dek5dEIsGoUaOQlpamX6fRaJCWlobIyMhO7UOtViMzMxPe3r3zi114tw5X7igg4rTz9BHSWVNaSsfbFbW4Ie99rY68l40rV67Ev//9b+zcuRPXr1/Ha6+9htraWixcuBAAsGDBAiQmJuq3X7t2LX788Ufcvn0bFy5cwEsvvYSCggIsXryYrx/BrA5mac+6nhzUD272Up6jIUJiL7VCdEvp2BvvdbTiO4AXXngBFRUVePvttyGXyzFixAikpKToL+IXFhZCJHqQY+/fv49XX30VcrkcLi4uGDVqFE6fPo2QkBC+fgSz0g9/QyUjMUJshDd+vFaG5CulWDltCDiO4zskk+FYbyyGH6O6uhpOTk5QKBT8Xf+qrQXs7bX/VyoBO7s2Nyu4W4uJ/zwOEQdk/HUq+tGZVyu1jbWwT9K+l8pEJewkbb+XfZVS1YxR76ZC1azBwdefRohPO5/5Tn4me0Jnv6O8l42kfboL9VGD3ShxEaPYS60wKdgDQO+7XYiSlwWjkpGYQuxDk3P0pkKLkpeFyqusxdWSaohFHGJCqZWRGG/yUA/IrEXIv1uHqyXCvMOkLZS8LNRBfcnYD652Ep6jIUJmJ7XC5KG60rH3tDpS8rJQB2iSDWJCuvkOelPpSMnLAt2uUOJ6KZWMxHR0pWPhvTpkFfeO0pH3fl6kNV3JOD7QDS49VDIyxtDc3Ay1Wt0jxzMVVaMK/nb+2v83qCDWiHmOqGvEYjGsrKzM3v/KVmKFKUM9kZxZiuTMUoT7Opn1eD2BkpcF0pWMs3qoZGxsbERpaSnq6up65HimpGEabBm/BQBQeqcUIk54xYStrS28vb0hkZj3D1VshHdL8irBm88EC77DKiUvC5NbrsQNeQ2sRBymh5p/kEWNRoO8vDyIxWL4+PhAIpEI6kOt1qhRX1kPAAhwC4BYJJwzL8YYGhsbUVFRgby8PAQFBRncTWJqk4I9YGMtRtG9emQWKxDh62y2Y/UESl4WRlcyPhXkBmdb85eMjY2N0Gg08PPzg62trdmPZ2pqjVr/KZbJZIJKXgBgY2MDa2trFBQUoLGx0azzMdhIxJgyzAMHrpQi+Uqp4JOX8M6xe7lknloZzfkXnzxeT773+mFyekGrI31iLcjNshpkl9XAWsxhegi1MhLTiw72gK1EjOKqely+o+A7nG6h5GVBdB0Inw5yh5OtNc/RkN5IZi3G1GHaa6lCH2GVkpcF0ZWMM6ljKjGj3tJhlZKXhcgpq8HNciWsxRymhfT+qdxMRS6X459/+ydmR86GrY0t/Pz88OyzzxqMzttd0dHRWLFihcn2x7foYHfYScQoUTTgYlEV3+EYjZKXhdCddU0IcoeTDZWMnZGfn4+xY8Yi41QGlv9tOS5dvoSUlBRMmjQJCQkJfIfXSmNjI98hAGgpHVv+QB4U8LyOlLwsAGPMoibZYIyhtrG2x5euljC///3vwXEcdibvxOTYyRgyZAhCQ0OxcuVK/PLLLwCAqqoqLF68GO7u7nB0dMTkyZNx+fJl/T5Wr16NESNG4D//+Q8CAgLg5OSEF198ETU12jHff/vb3+LEiRP4+OOPwXEcOI5Dfn4+ACArKwszZsyAvb09PD098fLLL6OyslK/7+joaCxbtgwrVqyAm5sbYmJiuvmbMR1da/bBzFJoNMIsHamflwXIKVMit1wJiVik/4vIp7qmOv3opD2pKyOh3rt3DykpKXj3f96Fja1Nq+ednZ0BAHPmzIGNjQ0OHToEJycnbN26FVOmTEFOTg5cXV0BALdu3cL+/ftx4MAB3L9/H3PnzsU//vEPrFu3Dh9//DFycnIQFhaGtWvXAgDc3d1RVVWFyZMnY/HixdiwYQPq6+vx5ptvYu7cuTh69Kg+jp07d+K1117DqVOnuvnumNaEIe6wl1rpS8dRbsIbuYSSlwXQtfpMGOIORxmVjJ2Rm5sLxhiGBg9td5uTJ0/i7NmzKC8vh1SqHYn2f//3f7F//37s3bsXS5YsAaC9y2DHjh1wcHAAALz88stIS0vDunXr4OTkBIlEAltbW4Pp+DZt2oSRI0fivffe06/77LPP4Ofnh5ycHAwZMgQAEBQUhA8++MDkP393yazFmBbiiX0Xi5F8pRSjJvvzHVKXUfLiGWMMB1pKxlkWUDICgK21LZSJPT/Lsq1153v4d6bEvHz5MpRKJfr162ewvr6+Hrdu3dI/DggI0CcuAPD29kZ5eXmH+z527Bjs7Vufod66dUufvEaNGtVhnHyJDffGvovFOJhZir9FDxDcNSRKXjzLLqvB7YpaSKxEmDLMg+9wAAAcx1n8RBZBQUHgOA43sm/AP7LtswalUglvb28cP3681XO6shLQTlz8MI7joNFoHnt8pVKJZ599Fu+//36r5x6eQ9SOx4ksOvL0EDc4SK0gr27ApaL7eILvgLqIkhfPDmfJAQDRQ9zhQCVjp7m6uiImJgab/7UZE+ZOaHXdq6qqCk888QTkcjmsrKwQEBBg9LEkEkmroYKeeOIJfPPNNwgICICVlTC/RlIrMaaFeuLbC8U4lCUXXPIS2plir5NyVZu8LKGVUWg++eQTqNVqxMfG42jyUdy8eRPXr1/H//3f/yEyMhJTp05FZGQk4uLi8OOPPyI/Px+nT5/GX//6V5w7d67TxwkICMCZM2eQn5+PyspKaDQaJCQk4N69e5g3bx4yMjJw69YtHD58GAsXLhTUmGi6VsfDLZ9DIaHkxbP8yrqWkpH/VkahGTRoEDLOZWB01Gh8tPYjDI8YjmnTpiEtLQ2bN28Gx3E4ePAgJkyYgIULF2LIkCF48cUXUVBQoJ/UuDP+9Kc/QSwWIyQkBO7u7igsLISPjw9OnToFtVqN6dOnIzw8HCtWrICzs7OgbnJ/KsgNDjIrVNRYRh+0LmEWYNOmTczf359JpVI2duxYdubMmcdu/9VXX7Hg4GAmlUpZWFgYS05O7vSxFAoFA8AUCkV3wzaeUskYwBjAhr6xly35IoO3UOrr69m1a9dYfX09bzF0R7O6mWUUZ7CM4gzWrG7mOxyj8P07WLnnEhv6xl79Z5IplbzEodPZ7yjvfyL27NmDlStX4p133sGFCxcwfPhwxMTEtNvac/r0acybNw+LFi3CxYsXERcXh7i4OGRlZZk8towM4MMPtf+aS3ONFLERPuY7ACEdeLSVu7jY9Me4cwc4dkz7r8n0UDJt19ixY1lCQoL+sVqtZj4+PiwpKanN7efOnctiY2MN1o0bN44tXbq0U8frbFaPj3/whwjQPjaVHZ88OPOyRQ37ZDN/Zwx8/9XvLjrz6j5Vk5oFzjij/0zac0q2fbvp9r99O2MikXb3IhHrcN+d/Y7y2kzS2NiI8+fPIzExUb9OJBJh6tSpSE9Pb/M16enpWLlypcG6mJgY7N+/v83tVSoVVCqV/nF1dcczp2RkADt3MgAPhkPeuZMh2/UinAZ0r/9TQ5UEGWvCEK9fw+H1ZWI8Nwvw9e3WrgkxSrlchNJDIfrHGgYsfpXhP0VnIHPu3rWwhioJTqwdBzDtd0mjAZYuBWJiuv955zV5VVZWQq1Wt7p46unpiRs3brT5Grlc3ub2cnnbrSVJSUlYs2ZNl+L6+Wfg4cSlxSH7shSOku7dyNpQ0K/VvtVqIDeXkhfhx82bAHv088443MhhkA2o6da+Gwr66ROXjqk+78LsoNIFiYmJBmdq1dXV8PPze+xrnn4aAAzPvACG93/vg2ER3WsVLJeLEL+HaXffQiwGAgO7tdtuYwIe10no+H7vg4IAEQeDz6RIxPCv14Lh4fX4zrodKZeLMOcrBo3mwXfJVJ93XpOXm5sbxGIxysrKDNaXlZUZ3Ef2MC8vry5tL5VK9fe1ddaYMUB8PIedOx+si4/n8OqvXbq0nzYFATWbaoGWEVvEImDrVv7OunS9y+vq6mBj0/oGZ2J+uinnHu3p31N8fYFNm2Dwmdy2jcPzT7t2f+dBwLZt2lJRrdYmLlN93jnGc9ofN24cxo4di40bNwLQ3iQ7YMAALFu2DH/5y19abf/CCy+grq4OP/zwg35dVFQUIiIisGXLlg6PV11dDScnJygUCjg6Oj5224wM4NQpYPx4bUIzmdpaoOWeuOJsJfoP4fcWktLSUlRVVcHDwwO2traCm/rseuV1AMAwt2GCmj2IMYa6ujqUl5fD2dnZ4LaiHmfmz+SdO9pSMTCw48TV2e8o72XjypUrER8fj9GjR2Ps2LH46KOPUFtbi4ULFwIAFixYgP79+yMpKQkAsHz5ckycOBHr169HbGwsvvzyS5w7dw7btm0zeWxjxpg4abWhf3/z7r8zdGetHd2MbIk0TINKhXYMrXxlviAnnXV2dm63cuCDOT6Tvr6mry54T14vvPACKioq8Pbbb0Mul2PEiBFISUnRX5QvLCw06LEcFRWF3bt3429/+xveeustBAUFYf/+/QgLC+PrRxA8juPg7e0NDw8PNDU18R1Ol9Q11iH2YCwA4MKSC7CVCGvuSWtra4jFwjlbtCS8l409rStlo9k8dIoOpRKw4JEHLF1tY61+4MSuDGZIHmFBn8nOfkeFd45NCCGg5EUIEShKXoQQQeL9gn1P013i68xtQmZTW/vg/9XV2g4wxCi1jbVAg/b/1dXVUEvovTSKBX0mdd/Nji7H97kL9nfu3Omwhz0hhH9FRUXwfUz/ij6XvDQaDUpKSuDg4NBhZ0zdrURFRUX8tUwaiWLnB8XefYwx1NTUwMfH57EDO/a5slEkEj02m7fF0dFRcB9EHYqdHxR79zg5OXW4DV2wJ4QIEiUvQoggUfJ6DKlUinfeeafLo1JYAoqdHxR7z+lzF+wJIb0DnXkRQgSJkhchRJAoeRFCBImSFyFEkCh5teOTTz5BQEAAZDIZxo0bh7Nnz/IdUoeSkpIwZswYODg4wMPDA3FxccjOzuY7LKP84x//AMdxWLFiBd+hdEpxcTFeeukl9OvXDzY2NggPD8e5c+f4DqtDarUaf//73zFw4EDY2Nhg8ODBePfdd3mfFKQzKHm1oauzeFuKEydOICEhAb/88gtSU1PR1NSE6dOno/bhm24FICMjA1u3bkVERATfoXTK/fv3MX78eFhbW+PQoUO4du0a1q9fDxcXE0zYYmbvv/8+Nm/ejE2bNuH69et4//338cEHH+jnlLBoxs+D23t1dRZvS1VeXs4AsBMnTvAdSqfV1NSwoKAglpqayiZOnMiWL1/Od0gdevPNN9lTTz3FdxhGiY2NZa+88orBuueff57Nnz+fp4g6j868HqGbxXvq1Kn6dR3N4m2pFAoFAMDV1QRTWPWQhIQExMbGGrz/lu7777/H6NGjMWfOHHh4eGDkyJH497//zXdYnRIVFYW0tDTk5OQAAC5fvoyTJ09ixowZPEfWsT53Y3ZHjJnF2xJpNBqsWLEC48ePF8zkJF9++SUuXLiAjIwMvkPpktu3b2Pz5s1YuXIl3nrrLWRkZOD111+HRCJBfHw83+E91l/+8hdUV1dj6NChEIvFUKvVWLduHebPn893aB2i5NVLJSQkICsrCydPnuQ7lE4pKirC8uXLkZqaCplMxnc4XaLRaDB69Gi89957AICRI0ciKysLW7Zssfjk9dVXX+G///0vdu/ejdDQUFy6dAkrVqyAj4+PxcdO17weoVKpmFgsZvv27TNYv2DBAvbcc8/xE1QXJSQkMF9fX3b79m2+Q+m0ffv2MQBMLBbrFwCM4zgmFotZc3Mz3yG2a8CAAWzRokUG6/71r38xHx8fniLqPF9fX7Zp0yaDde+++y4LDg7mKaLOo2tej5BIJBg1ahTS0tL06zQaDdLS0hAZGcljZB1jjGHZsmXYt28fjh49ioEDB/IdUqdNmTIFmZmZuHTpkn4ZPXo05s+fj0uXLln03Ibjx49v1SUlJycH/v7+PEXUeXV1da0G/BOLxdBoNDxF1AV8Z09L9OWXXzKpVMp27NjBrl27xpYsWcKcnZ2ZXC7nO7THeu2115iTkxM7fvw4Ky0t1S91dXV8h2YUobQ2nj17lllZWbF169axmzdvsv/+97/M1taW7dq1i+/QOhQfH8/69+/PDhw4wPLy8ti3337L3Nzc2KpVq/gOrUOUvNqxceNGNmDAACaRSNjYsWPZL7/8wndIHQLQ5vL555/zHZpRhJK8GGPshx9+YGFhYUwqlbKhQ4eybdu28R1Sp1RXV7Ply5ezAQMGMJlMxgYNGsT++te/MpVKxXdoHaIhcQghgkTXvAghgkTJixAiSJS8CCGCRMmLECJIlLwIIYJEyYsQIkiUvAghgkTJixAiSJS8CCGCRMmL9CoFBQWwsbGBUqnkOxRiZpS8SK/y3XffYdKkSbC3t+c7FGJmlLyIRYqOjsYf/vAHrFixAi4uLvD09MS///1v1NbWYuHChXBwcEBgYCAOHTpk8LrvvvsOzz33HACA47hWS0BAAA8/DTEHSl7EYu3cuRNubm44e/Ys/vCHP+C1117DnDlzEBUVhQsXLmD69Ol4+eWXUVdXBwCoqqrCyZMn9cmrtLRUv+Tm5iIwMBATJkzg80ciJkSjShCLFB0dDbVajZ9//hmAdn5BJycnPP/88/jiiy8AAHK5HN7e3khPT8eTTz6J3bt3Y8OGDa3GwGeM4de//jUKCwvx888/w8bGpsd/HmJ6NIY9sVgPz9soFovRr18/hIeH69fpJknRzaf5cMn4sLfeegvp6ek4d+4cJa5ehMpGYrGsra0NHnMcZ7CO4zgA2mG6GxsbkZKS0ip57dq1Cxs2bMC+ffvQv39/8wdNegwlL9IrHD9+HC4uLhg+fLh+XXp6OhYvXoytW7fiySef5DE6Yg5UNpJe4fvvvzc465LL5fjVr36FF198ETExMZDL5QC05ae7uztfYRITojMv0is8mrxu3LiBsrIy7Ny5E97e3vplzJgxPEZJTIlaG4ngXbhwAZMnT0ZFRUWr62Sk96IzLyJ4zc3N2LhxIyWuPobOvAghgkRnXoQQQaLkRQgRJEpehBBBouRFCBEkSl6EEEGi5EUIESRKXoQQQaLkRQgRJEpehBBB+v9q8P+oGzNbNgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def plot_centroid(peak, mz_array, int_array):\n",
    "\n",
    "    plt.figure(figsize=(3, 3))\n",
    "\n",
    "    start, center, end = peak\n",
    "\n",
    "    centroid = get_centroid((start, center, end), mz_array, int_array)\n",
    "\n",
    "    plt.plot(mz_array[start - 2 : end + 2], int_array[start - 2 : end + 2])\n",
    "    plt.plot(mz_array[start - 2 : end + 2], int_array[start - 2 : end + 2], \"b.\")\n",
    "    plt.axvline(mz_array[start], color=\"r\")\n",
    "    plt.axvline(mz_array[end], color=\"r\")\n",
    "\n",
    "    plt.axvline(centroid[0], color=\"g\", label='Center')\n",
    "    plt.title(\"Centroid\")\n",
    "    plt.xlabel(\"m/z\")\n",
    "    plt.ylabel(\"Intensity\")\n",
    "    plt.legend()\n",
    "\n",
    "    plt.show()\n",
    "    \n",
    "int_array = np.array([0,0,0,1,2,3,3,2,1,0,0,1,2,3,2,1,0])\n",
    "mz_array = np.array([0,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])\n",
    "\n",
    "peak = (2, 6, 9)\n",
    "    \n",
    "plot_centroid(peak, mz_array, int_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "@njit\n",
    "def centroid_data(\n",
    "    mz_array: np.ndarray,\n",
    "    int_array: np.ndarray\n",
    ") -> tuple:\n",
    "    \"\"\"Estimate centroids and intensities from profile data.\n",
    "\n",
    "    Args:\n",
    "        mz_array (np.ndarray): An array with mz values.\n",
    "        int_array (np.ndarray): An array with intensity values.\n",
    "\n",
    "    Returns:\n",
    "        tuple: A tuple of the form (mz_array_centroided, int_array_centroided)\n",
    "    \"\"\"\n",
    "    peaks = get_peaks(int_array)\n",
    "\n",
    "    mz_array_centroided = np.zeros(len(peaks))\n",
    "    int_array_centroided = np.zeros(len(peaks))\n",
    "\n",
    "\n",
    "    for i in range(len(peaks)):\n",
    "        mz_, int_ = get_centroid(peaks[i], mz_array, int_array)\n",
    "        mz_array_centroided[i] = mz_\n",
    "        int_array_centroided[i] = int_\n",
    "\n",
    "    return mz_array_centroided, int_array_centroided"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trimming spectra to retain the `n` most intense peaks\n",
    "\n",
    "`get_most_abundant`: In order to save spectra in a more memory-efficient form, we only keep the n most abundant peaks. This allows us to save data in a fast, accessible matrix format. \n",
    "\n",
    "`get_local_intensity`: This calculates the local intensity to get local maxima."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "from alphapept.chem import calculate_mass\n",
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "from numba.typed import List\n",
    "from numba import njit\n",
    "import gzip\n",
    "import sys\n",
    "import os\n",
    "import logging\n",
    "\n",
    "@njit\n",
    "def get_local_intensity(intensity, window=10):\n",
    "    \"\"\"\n",
    "    Calculate the local intensity for a spectrum.\n",
    "\n",
    "    Args:\n",
    "        intensity (np.ndarray): An array with intensity values.\n",
    "        window (int): Window Size\n",
    "    Returns:\n",
    "        nop.ndarray: local intensity\n",
    "    \"\"\"\n",
    "\n",
    "    local_intensity = np.zeros(len(intensity))\n",
    "\n",
    "    for i in range(len(intensity)):\n",
    "        start = max(0, i-window) \n",
    "        end = min(len(intensity), i+window)\n",
    "        local_intensity[i] = intensity[i]/np.max(intensity[start:end])\n",
    "        \n",
    "    return local_intensity\n",
    "\n",
    "def get_most_abundant(\n",
    "    mass: np.ndarray,\n",
    "    intensity: np.ndarray,\n",
    "    n_max: int,\n",
    "    window: int = 10,\n",
    ") -> tuple:\n",
    "    \"\"\"Returns the n_max most abundant peaks of a spectrum.\n",
    "\n",
    "    Args:\n",
    "        mass (np.ndarray): An array with mz values.\n",
    "        intensity (np.ndarray): An array with intensity values.\n",
    "        n_max (int): The maximum number of peaks to retain.\n",
    "            Setting `n_max` to -1 returns all peaks.\n",
    "        window (int): Use local maximum in a window\n",
    "\n",
    "    Returns:\n",
    "        tuple: the filtered mass and intensity arrays.\n",
    "\n",
    "    \"\"\"\n",
    "    if n_max == -1:\n",
    "        return mass, intensity\n",
    "    if len(mass) < n_max:\n",
    "        return mass, intensity\n",
    "    else:\n",
    "        \n",
    "        if window > 0:\n",
    "            sortindex = np.argsort(get_local_intensity(intensity, window))[::-1][:n_max]\n",
    "        else:\n",
    "            sortindex = np.argsort(intensity)[::-1][:n_max]\n",
    "        \n",
    "        sortindex.sort()\n",
    "\n",
    "    return mass[sortindex], intensity[sortindex]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading raw input\n",
    "\n",
    "Input can be read from three different files formats:\n",
    "\n",
    "* `.raw` files describing `Thermo` raw data.\n",
    "* `.d` folders describing `Bruker` raw data.\n",
    "* `.mzML` files describing generic input.\n",
    "* `.wiff` files for `SciEx` raw data.\n",
    "\n",
    "All reading functions return a `query_data` dictionary similar to the following structure:\n",
    "\n",
    "```\n",
    "query_data = {\n",
    "    \"scan_list_ms1\":   np.array(...)\n",
    "    \"rt_list_ms1\":     np.array(...)\n",
    "    \"mass_list_ms1\":   np.array(...)\n",
    "    \"int_list_ms1\":    np.array(...)\n",
    "    \"ms_list_ms1\":     np.array(...)\n",
    "    \"scan_list_ms2\":   np.array(...)\n",
    "    \"rt_list_ms2\":     np.array(...)\n",
    "    \"mass_list_ms2\":   mass_list_ms2\n",
    "    \"int_list_ms2\":    int_list_ms2\n",
    "    \"ms_list_ms2\":     np.array(...)\n",
    "    \"prec_mass_list2\": np.array(...)\n",
    "    \"mono_mzs2\":       np.array(...)\n",
    "    \"charge2\":         np.array(...)\n",
    "}\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thermo\n",
    "\n",
    "To read Thermo files, AlphaPept uses the `pyrawfilereader` package, a Python implementation of the commonly used `rawfilereader` tool. By using the custom python version, Thermo files can be read without having to install `MSFileReader`.\n",
    "\n",
    "The user can pass an additional flag `use_profile_ms1`. This will then use the profile data which is not centroided already an peform centroiding. Note that this will lead to slightly different intensities, as the centroided data uses the apex and the centroid algorithm the summed intensity.\n",
    "\n",
    "The import of files is done now via `alpharaw`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def load_thermo_raw(\n",
    "    raw_file_name: str,\n",
    "    n_most_abundant: int,\n",
    "    use_profile_ms1: bool = False,\n",
    "    callback: callable = None,\n",
    ") -> tuple:\n",
    "    \"\"\"Load raw thermo data as a dictionary.\n",
    "\n",
    "    Args:\n",
    "        raw_file_name (str): The name of a Thermo .raw file.\n",
    "        n_most_abundant (int): The maximum number of peaks to retain per MS2 spectrum.\n",
    "        use_profile_ms1 (bool): Use profile data or centroid it beforehand. Defaults to False.\n",
    "        callback (callable): A function that accepts a float between 0 and 1 as progress. Defaults to None.\n",
    "\n",
    "    Returns:\n",
    "        tuple: A dictionary with all the raw data and a string with the acquisition_date_time\n",
    "\n",
    "    \"\"\"\n",
    "    from alphapept.pyrawfilereader import RawFileReader\n",
    "    rawfile = RawFileReader(raw_file_name)\n",
    "\n",
    "    spec_indices = range(rawfile.FirstSpectrumNumber, rawfile.LastSpectrumNumber + 1)\n",
    "\n",
    "    scan_list = []\n",
    "    rt_list = []\n",
    "    mass_list = []\n",
    "    int_list = []\n",
    "    ms_list = []\n",
    "    prec_mzs_list = []\n",
    "    mono_mzs_list = []\n",
    "    charge_list = []\n",
    "\n",
    "    for idx, i in enumerate(spec_indices):\n",
    "        try:\n",
    "            ms_order = rawfile.GetMSOrderForScanNum(i)\n",
    "            rt = rawfile.RTFromScanNum(i)\n",
    "\n",
    "            if ms_order == 2:\n",
    "                prec_mz = rawfile.GetPrecursorMassForScanNum(i, 0)\n",
    "\n",
    "                mono_mz, charge = rawfile.GetMS2MonoMzAndChargeFromScanNum(i)\n",
    "            else:\n",
    "                prec_mz, mono_mz, charge = 0,0,0\n",
    "\n",
    "            if use_profile_ms1:\n",
    "                if ms_order == 2:\n",
    "                    masses, intensity = rawfile.GetCentroidMassListFromScanNum(i)\n",
    "                    masses, intensity = get_most_abundant(masses, intensity, n_most_abundant)\n",
    "                else:\n",
    "                    masses, intensity = rawfile.GetProfileMassListFromScanNum(i)\n",
    "                    masses, intensity = centroid_data(masses, intensity)\n",
    "\n",
    "            else:\n",
    "                masses, intensity = rawfile.GetCentroidMassListFromScanNum(i)\n",
    "                if ms_order == 2:\n",
    "                    masses, intensity = get_most_abundant(masses, intensity, n_most_abundant)\n",
    "\n",
    "            scan_list.append(i)\n",
    "            rt_list.append(rt)\n",
    "            mass_list.append(np.array(masses))\n",
    "            int_list.append(np.array(intensity, dtype=np.int64))\n",
    "            ms_list.append(ms_order)\n",
    "            prec_mzs_list.append(prec_mz)\n",
    "            mono_mzs_list.append(mono_mz)\n",
    "            charge_list.append(charge)\n",
    "        except KeyboardInterrupt as e:\n",
    "            raise e\n",
    "        except SystemExit as e:\n",
    "            raise e\n",
    "        except Exception as e:\n",
    "            logging.info(f\"Bad scan={i} in raw file '{raw_file_name}': {e}\")\n",
    "\n",
    "        if callback:\n",
    "            callback((idx+1)/len(spec_indices))\n",
    "\n",
    "    scan_list_ms1 = [scan_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    rt_list_ms1 = [rt_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    mass_list_ms1 = [mass_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    int_list_ms1 = [int_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    ms_list_ms1 = [ms_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "\n",
    "    scan_list_ms2 = [scan_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    rt_list_ms2 = [rt_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    mass_list_ms2 = [mass_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    int_list_ms2 = [int_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    ms_list_ms2 = [ms_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    mono_mzs2 = [mono_mzs_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    charge2 = [charge_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "\n",
    "    prec_mass_list2 = [\n",
    "        calculate_mass(mono_mzs_list[i], charge_list[i])\n",
    "        for i, _ in enumerate(ms_list)\n",
    "        if _ == 2\n",
    "    ]\n",
    "\n",
    "    check_sanity(mass_list)\n",
    "\n",
    "    query_data = {}\n",
    "\n",
    "    query_data[\"scan_list_ms1\"] = np.array(scan_list_ms1)\n",
    "    query_data[\"rt_list_ms1\"] = np.array(rt_list_ms1)\n",
    "    query_data[\"mass_list_ms1\"] = np.array(mass_list_ms1, dtype=object)\n",
    "    query_data[\"int_list_ms1\"] = np.array(int_list_ms1, dtype=object)\n",
    "    query_data[\"ms_list_ms1\"] = np.array(ms_list_ms1)\n",
    "\n",
    "    query_data[\"scan_list_ms2\"] = np.array(scan_list_ms2)\n",
    "    query_data[\"rt_list_ms2\"] = np.array(rt_list_ms2)\n",
    "    query_data[\"mass_list_ms2\"] = mass_list_ms2\n",
    "    query_data[\"int_list_ms2\"] = int_list_ms2\n",
    "    query_data[\"ms_list_ms2\"] = np.array(ms_list_ms2)\n",
    "    query_data[\"prec_mass_list2\"] = np.array(prec_mass_list2)\n",
    "    query_data[\"mono_mzs2\"] = np.array(mono_mzs2)\n",
    "#     TODO: Refactor charge2 to be consistent: charge_ms2\n",
    "    query_data[\"charge2\"] = np.array(charge2)\n",
    "\n",
    "    acquisition_date_time = rawfile.GetCreationDate()\n",
    "\n",
    "    rawfile.Close()\n",
    "\n",
    "    return query_data, acquisition_date_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bruker\n",
    "\n",
    "To access Bruker files, AlphaPept relies on the external `timsdata` library from Bruker (available in the `alphatims\\ext` folder, licenses are applicable). Unfortunately, these libraries are only available on Windows and Linux. As a result, the reading of raw data is not available on macOS. However, once raw data is converted to `.ms_data.hdf` output, other workflow steps (besides feature feating) are possible without problems on macOS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def load_bruker_raw(raw_file, most_abundant, callback=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Load bruker raw file and extract spectra\n",
    "    \"\"\"\n",
    "    import alphatims.bruker\n",
    "    from .constants import mass_dict\n",
    "    from .io import list_to_numpy_f32, get_most_abundant\n",
    "\n",
    "    data = alphatims.bruker.TimsTOF(raw_file)\n",
    "    prec_data = data.precursors\n",
    "    frame_data = data.frames\n",
    "    frame_data = frame_data.set_index('Id')\n",
    "    \n",
    "    import sqlalchemy as db\n",
    "    import pandas as pd\n",
    "\n",
    "    tdf = os.path.join(raw_file, 'analysis.tdf')\n",
    "    engine = db.create_engine('sqlite:///{}'.format(tdf))\n",
    "\n",
    "    global_metadata = pd.read_sql_table('GlobalMetadata', engine)\n",
    "    global_metadata = global_metadata.set_index('Key').to_dict()['Value']\n",
    "    acquisition_date_time = global_metadata['AcquisitionDateTime']\n",
    "    \n",
    "\n",
    "    M_PROTON = mass_dict['Proton']\n",
    "\n",
    "    prec_data['Mass'] = prec_data['MonoisotopicMz'].values * prec_data['Charge'].values - prec_data['Charge'].values*M_PROTON\n",
    "\n",
    "    query_data = {}\n",
    "\n",
    "    query_data['prec_mass_list2'] = prec_data['Mass'].values\n",
    "    query_data['prec_id2'] = prec_data['Id'].values\n",
    "    query_data['mono_mzs2'] = prec_data['MonoisotopicMz'].values\n",
    "    query_data['rt_list_ms2'] = frame_data.loc[prec_data['Parent'].values]['Time'].values / 60 #convert to minutes\n",
    "    query_data['scan_list_ms2'] = prec_data['Parent'].values\n",
    "    query_data['charge2'] = prec_data['Charge'].values\n",
    "\n",
    "    query_data['mobility2'] = data.mobility_values[\n",
    "        data.precursors.ScanNumber.values.astype(np.int64)\n",
    "    ]\n",
    "    (\n",
    "        spectrum_indptr,\n",
    "        spectrum_tof_indices,\n",
    "        spectrum_intensity_values,\n",
    "    ) = data.index_precursors(\n",
    "        centroiding_window=5,\n",
    "        keep_n_most_abundant_peaks=most_abundant\n",
    "    )\n",
    "    # TODO: Centroid spectra and trim\n",
    "    query_data[\"alphatims_spectrum_indptr_ms2\"] = spectrum_indptr[1:]\n",
    "    query_data[\"alphatims_spectrum_mz_values_ms2\"] = data.mz_values[spectrum_tof_indices]\n",
    "    query_data[\"alphatims_spectrum_intensity_values_ms2\"] = spectrum_intensity_values\n",
    "\n",
    "    return query_data, acquisition_date_time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For `ccs` (i.e., ion mobility) values, we need additional functions from the Bruker library. As the live feature-finder might not be able to determine some charge values, it is intended to perform this calculation at a later stage once we have charge values from the post-processing feature finder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "import alphapept\n",
    "import numpy as np\n",
    "\n",
    "def one_over_k0_to_CCS(\n",
    "    one_over_k0s: np.ndarray,\n",
    "    charges: np.ndarray,\n",
    "    mzs: np.ndarray,\n",
    ") -> np.ndarray:\n",
    "    \"\"\"Retrieve collisional cross section (CCS) values from (mobility, charge, mz) arrays.\n",
    "\n",
    "    Args:\n",
    "        one_over_k0s (np.ndarray): The ion mobilities (1D-np.float).\n",
    "        charges (np.ndarray): The charges (1D-np.int).\n",
    "        mzs (np.ndarray): The mz values (1D-np.float).\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: The CCS values.\n",
    "\n",
    "    \"\"\"\n",
    "    from alphapept.ext.bruker import timsdata\n",
    "\n",
    "    ccs = np.empty(len(one_over_k0s))\n",
    "    ccs[:] = np.nan\n",
    "\n",
    "    for idx, (one_over, charge, mz) in enumerate(zip(one_over_k0s, charges, mzs)):\n",
    "        try:\n",
    "            ccs[idx] = timsdata.oneOverK0ToCCSforMz(one_over, int(charge), mz)\n",
    "        except ValueError:\n",
    "            pass\n",
    "    return ccs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to availbility of Bruker libraries, this can only be tested on Windows and Linux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "def test_one_over_k0_to_CCS():\n",
    "    one_over_k0_to_CCS([1], [1], [1])\n",
    "\n",
    "# NOTE, test was not enable for a long time and might be outdated by now\n",
    "# import sys\n",
    "# if sys.platform[:6] != \"darwin\":\n",
    "#     test_one_over_k0_to_CCS()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  `.wiff` files\n",
    "\n",
    "Similar to Thermo, we now use alpharaw to read `.wiff`-files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "\n",
    "from alpharaw.wrappers.alphapept_wrapper import parse_msdata_to_alphapept\n",
    "from alphapept.chem import calculate_mass\n",
    "\n",
    "def import_sciex_as_alphapept(sciex_wiff):\n",
    "    from alpharaw.sciex import SciexWiffData\n",
    "    raw_data = SciexWiffData(centroided=False)\n",
    "    raw_data.import_raw(sciex_wiff)\n",
    "    return parse_msdata_to_alphapept(raw_data)\n",
    "\n",
    "def load_sciex_raw(\n",
    "    raw_file_name: str,\n",
    "    n_most_abundant: int,\n",
    "    callback: callable = None,\n",
    ") -> tuple:\n",
    "    \"\"\"Load Sciex data as a dictionary.\n",
    "\n",
    "    Args:\n",
    "        raw_file_name (str): The name of a Sciex .wiff file.\n",
    "        n_most_abundant (int): The maximum number of peaks to retain per MS2 spectrum.\n",
    "        use_profile_ms1 (bool): Use profile data or centroid it beforehand. Defaults to False.\n",
    "        callback (callable): A function that accepts a float between 0 and 1 as progress. Defaults to None.\n",
    "\n",
    "    Returns:\n",
    "        tuple: A dictionary with all the raw data and a string with the acquisition_date_time\n",
    "\n",
    "    Note: SciEx file have one main file but multiple additional files which need to be present in the same location.\n",
    "    \"\"\"\n",
    "    query_data, acquisition_date_time = import_sciex_as_alphapept(raw_file_name)\n",
    "\n",
    "    query_data['prec_mass_list2'] = np.array([\n",
    "        calculate_mass(query_data['mono_mzs2'][i],  query_data['charge2'][i])\n",
    "        for i in range(len(query_data['mono_mzs2']))])\n",
    "\n",
    "    # Centroid data\n",
    "    # ToDo: Should we insert zero intensity peaks?\n",
    "    for index in ['1','2']:\n",
    "        masses_ = []\n",
    "        intensities_ = []\n",
    "        for j in range(len(query_data[f'mass_list_ms{index}'])):\n",
    "            mass = np.array(query_data[f'mass_list_ms{index}'][j])\n",
    "            ints = np.array(query_data[f'int_list_ms{index}'][j])\n",
    "            sortindex = np.argsort(mass)\n",
    "            mass = mass[sortindex]\n",
    "            ints = ints[sortindex]\n",
    "            mass_, ints_ = centroid_data(np.array(mass), np.array(ints))\n",
    "            if index == '2':\n",
    "                mass_, ints_ = get_most_abundant(mass_, ints_, n_most_abundant)\n",
    "            masses_.append(mass_)\n",
    "            intensities_.append(ints_)\n",
    "\n",
    "        query_data[f'mass_list_ms{index}'] = np.array(masses_, dtype='object')\n",
    "        query_data[f'int_list_ms{index}'] = np.array(intensities_, dtype='object')\n",
    "\n",
    "    return query_data, acquisition_date_time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generic `.mzML` files\n",
    "\n",
    "To access `.mzML` files, we rely on the `pyteomics` package. For using an mzml format for performing a search, Peak Picking (data centroiding) should be applied to all MS levels of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def check_sanity(mass_list: np.ndarray) -> None:\n",
    "    \"\"\"Sanity check for mass list to make sure the masses are sorted.\n",
    "\n",
    "    Args:\n",
    "        mass_list (np.ndarray): The mz values (1D-np.float).\n",
    "\n",
    "    Raises:\n",
    "        ValueError: When the mz values are not sorted.\n",
    "\n",
    "    \"\"\"\n",
    "    if not all(\n",
    "        mass_list[0][i] <= mass_list[0][i + 1] for i in range(len(mass_list[0]) - 1)\n",
    "    ):\n",
    "        raise ValueError(\"Masses are not sorted.\")\n",
    "\n",
    "\n",
    "def extract_mzml_info(input_dict: dict) -> tuple:\n",
    "    \"\"\"Extract basic MS coordinate arrays from a dictionary.\n",
    "\n",
    "    Args:\n",
    "        input_dict (dict): A dictionary obtained by iterating over a Pyteomics mzml.read function.\n",
    "\n",
    "    Returns:\n",
    "        tuple: The rt, masses, intensities, ms_order, prec_mass, mono_mz, charge arrays retrieved from the input_dict.\n",
    "            If the `ms level` in the input dict does not equal 2, the charge, mono_mz and prec_mass will be equal to 0.\n",
    "\n",
    "    \"\"\"\n",
    "    from alphapept.chem import calculate_mass\n",
    "    rt = float(input_dict.get('scanList').get('scan')[0].get('scan start time'))  # rt_list_ms1/2\n",
    "    masses = input_dict.get('m/z array')\n",
    "    intensities = input_dict.get('intensity array').astype(int)\n",
    "    ms_order = input_dict.get('ms level')\n",
    "    prec_mass = mono_mz = charge = 0\n",
    "    if ms_order == 2:\n",
    "        try:\n",
    "            charge = int(input_dict.get('precursorList').get('precursor')[0].get('selectedIonList').get('selectedIon')[0].get(\n",
    "                'charge state'))\n",
    "        except TypeError:\n",
    "            charge = 0\n",
    "        mono_mz = input_dict.get('precursorList').get('precursor')[0].get('selectedIonList').get('selectedIon')[0].get(\n",
    "                'selected ion m/z')\n",
    "        prec_mass = calculate_mass(mono_mz, charge)\n",
    "    return rt, masses, intensities, ms_order, prec_mass, mono_mz, charge\n",
    "\n",
    "\n",
    "def load_mzml_data(\n",
    "    filename: str,\n",
    "    n_most_abundant: int,\n",
    "    callback: callable = None,\n",
    "    **kwargs\n",
    ") -> tuple:\n",
    "    \"\"\"Load data from an mzml file as a dictionary.\n",
    "\n",
    "    Args:\n",
    "        filename (str): The name of a .mzml file.\n",
    "        n_most_abundant (int): The maximum number of peaks to retain per MS2 spectrum.\n",
    "        callback (callable): A function that accepts a float between 0 and 1 as progress. Defaults to None.\n",
    "\n",
    "    Returns:\n",
    "        tuple: A dictionary with all the raw data, a string with the acquisition_date_time and a string with the vendor.\n",
    "\n",
    "    \"\"\"\n",
    "    from pyteomics import mzml\n",
    "    import os\n",
    "    import re\n",
    "    import logging\n",
    "    import datetime\n",
    "    import pathlib\n",
    "    np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)\n",
    "\n",
    "    try:\n",
    "        reader = mzml.read(filename, use_index=True)\n",
    "        spec_indices = np.array(range(1, len(reader) + 1))\n",
    "    except OSError:\n",
    "        logging.info('Could not open the file. Please, specify the correct path to the file.')\n",
    "        sys.exit(1)\n",
    "\n",
    "    scan_list = []\n",
    "    rt_list = []\n",
    "    mass_list = []\n",
    "    int_list = []\n",
    "    ms_list = []\n",
    "    prec_mzs_list = []\n",
    "    mono_mzs_list = []\n",
    "    charge_list = []\n",
    "    \n",
    "    vendor = \"Unknown\"\n",
    "    \n",
    "    for idx, i in enumerate(spec_indices):\n",
    "        try:\n",
    "            spec = next(reader)\n",
    "\n",
    "            if idx == 0:\n",
    "                ext = re.findall(r\"File:\\\".+\\.(\\w+)\\\"\", spec['spectrum title'])[0]\n",
    "                if ext.lower() == 'raw':\n",
    "                    vendor = \"Thermo\"\n",
    "\n",
    "            scan_list.append(i)\n",
    "            rt, masses, intensities, ms_order, prec_mass, mono_mz, charge = extract_mzml_info(spec)\n",
    "            \n",
    "            sortindex = np.argsort(masses)\n",
    "            \n",
    "            masses = masses[sortindex]\n",
    "            intensities = intensities[sortindex]\n",
    "            \n",
    "            if ms_order == 2:\n",
    "                masses, intensities = get_most_abundant(masses, intensities, n_most_abundant)\n",
    "            rt_list.append(rt)\n",
    "            \n",
    "            #Remove zero intensities\n",
    "            to_keep = intensities>0\n",
    "            masses = masses[to_keep]\n",
    "            intensities = intensities[to_keep]\n",
    "            \n",
    "            mass_list.append(masses)\n",
    "            int_list.append(intensities)\n",
    "            ms_list.append(ms_order)\n",
    "            prec_mzs_list.append(prec_mass)\n",
    "            mono_mzs_list.append(mono_mz)\n",
    "            charge_list.append(charge)\n",
    "            \n",
    "        except KeyboardInterrupt as e:\n",
    "            raise e\n",
    "        except SystemExit as e:\n",
    "            raise e\n",
    "        except Exception as e:\n",
    "            logging.info(f\"Bad scan={i} in mzML file '{filename}' {e}\")\n",
    "\n",
    "        if callback:\n",
    "            callback((idx+1)/len(spec_indices))\n",
    "\n",
    "    check_sanity(mass_list)\n",
    "\n",
    "    scan_list_ms1 = [scan_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    rt_list_ms1 = [rt_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    mass_list_ms1 = [mass_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    int_list_ms1 = [int_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "    ms_list_ms1 = [ms_list[i] for i, _ in enumerate(ms_list) if _ == 1]\n",
    "\n",
    "    scan_list_ms2 = [scan_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    rt_list_ms2 = [rt_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    mass_list_ms2 = [mass_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    int_list_ms2 = [int_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    ms_list_ms2 = [ms_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    prec_mass_list2 = [prec_mzs_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    mono_mzs2 = [mono_mzs_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "    charge_ms2 = [charge_list[i] for i, _ in enumerate(ms_list) if _ == 2]\n",
    "\n",
    "    prec_mass_list2 = [\n",
    "        calculate_mass(mono_mzs_list[i], charge_list[i])\n",
    "        for i, _ in enumerate(ms_list)\n",
    "        if _ == 2\n",
    "    ]\n",
    "\n",
    "    query_data = {}\n",
    "\n",
    "    query_data[\"scan_list_ms1\"] = np.array(scan_list_ms1)\n",
    "    query_data[\"rt_list_ms1\"] = np.array(rt_list_ms1)\n",
    "    query_data[\"mass_list_ms1\"] = np.array(mass_list_ms1)\n",
    "    query_data[\"int_list_ms1\"] = np.array(int_list_ms1)\n",
    "    query_data[\"ms_list_ms1\"] = np.array(ms_list_ms1)\n",
    "\n",
    "    query_data[\"scan_list_ms2\"] = np.array(scan_list_ms2)\n",
    "    query_data[\"rt_list_ms2\"] = np.array(rt_list_ms2)\n",
    "    query_data[\"mass_list_ms2\"] = mass_list_ms2\n",
    "    query_data[\"int_list_ms2\"] = int_list_ms2\n",
    "    query_data[\"ms_list_ms2\"] = np.array(ms_list_ms2)\n",
    "    query_data[\"prec_mass_list2\"] = np.array(prec_mass_list2)\n",
    "    query_data[\"mono_mzs2\"] = np.array(mono_mzs2)\n",
    "    query_data[\"charge2\"] = np.array(charge_ms2)\n",
    "\n",
    "    fname = pathlib.Path(filename)\n",
    "    acquisition_date_time = datetime.datetime.fromtimestamp(fname.stat().st_mtime).strftime('%Y-%m-%dT%H:%M:%S')\n",
    "\n",
    "    return query_data, acquisition_date_time, vendor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading other files\n",
    "\n",
    "Benchmarking proteomics software against each other is not straightforward as various naming conventions exist and different algorithms are implemented. In this section, we define some helper functions that allow us to read results from other tools and facilitate the comparison of different tools against AlphaPept."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading MaxQuant .xml settings file\n",
    "\n",
    "One of the most commonly used tools to analyze MS data is MaxQuant. AlphaPept reads MaxQuant .xml files as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "import xml.etree.ElementTree as ET\n",
    "\n",
    "def __extract_nested(child):\n",
    "    \"\"\"\n",
    "    Helper function to extract nested entries\n",
    "    \"\"\"\n",
    "    if len(child) > 0:\n",
    "        temp_dict = {}\n",
    "        for xx in child:\n",
    "            temp_dict[xx.tag] = __extract_nested(xx)\n",
    "        return temp_dict\n",
    "    else:\n",
    "        if child.text == 'True':\n",
    "            info = True\n",
    "        elif child.text == 'False':\n",
    "            info = False\n",
    "        else:\n",
    "            info = child.text\n",
    "        return info\n",
    "\n",
    "def extract_mq_settings(path: str) -> dict:\n",
    "    \"\"\"Function to return MaxQuant values as a dictionary for a given xml file.\n",
    "\n",
    "    Args:\n",
    "        path (str): File name of an xml file.\n",
    "\n",
    "    Returns:\n",
    "        dict: A dictionary with MaxQuant info.\n",
    "\n",
    "    Raises:\n",
    "        ValueError: When path is not a valid xml file.\n",
    "\n",
    "    \"\"\"\n",
    "    if not path.endswith('.xml'):\n",
    "        raise ValueError(\"Path {} is not a valid xml file.\".format(path))\n",
    "\n",
    "    tree = ET.parse(path)\n",
    "    root = tree.getroot()\n",
    "\n",
    "    mq_dict = {}\n",
    "\n",
    "    for child in root:\n",
    "\n",
    "        mq_dict[child.tag] = __extract_nested(child)\n",
    "\n",
    "    return mq_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'FastaFileInfo': {'fastaFilePath': 'testfile.fasta',\n",
       "  'identifierParseRule': '>([^\\\\s]*)',\n",
       "  'descriptionParseRule': '>(.*)',\n",
       "  'taxonomyParseRule': None,\n",
       "  'variationParseRule': None,\n",
       "  'modificationParseRule': None,\n",
       "  'taxonomyId': None}}"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mq_dict = extract_mq_settings('../testfiles/test_mqpar.xml')\n",
    "mq_dict['fastaFiles']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AlphaPept incorporates PTMs directly in amino acid sequences with lower case identifiers. Parsing MaxQuant sequences with PTMs is done with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def parse_mq_seq(peptide: str) -> str:\n",
    "    \"\"\"Replaces maxquant convention to alphapept convention.\n",
    "\n",
    "    ToDo: include more sequences\n",
    "\n",
    "    Args:\n",
    "        peptide (str): A peptide sequence from MaxQuant.\n",
    "\n",
    "    Returns:\n",
    "        str: A parsed peptide sequence compatible with AlphaPept.\n",
    "\n",
    "    \"\"\"\n",
    "    peptide = peptide[1:-1] #Remove _\n",
    "\n",
    "    peptide = peptide.replace('(Acetyl (Protein N-term))','a')\n",
    "    peptide = peptide.replace('M(Oxidation (M))','oxM')\n",
    "    peptide = peptide.replace('C','cC') #This is fixed and not indicated in MaxQuant\n",
    "\n",
    "    return peptide"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick test shows the results are correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'AFQPFFVELToxMPYSVIR'"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "parse_mq_seq('_AFQPFFVELTM(Oxidation (M))PYSVIR_')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For saving, we are currently relying on the hdf-container (see below).\n",
    "\n",
    "While we could, in principle, store the mz and int arrays as a list of variable length, this will come at a performance decrease. We, therefore, create an array of the dimensions of the n most abundant peaks and the number of spectra with the function `list_to_numpy_f32` and fill the unoccupied cells with `-1`. This allows an increase in accessing times at the cost of additional disk space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def list_to_numpy_f32(\n",
    "    long_list: list\n",
    ") -> np.ndarray:\n",
    "    \"\"\"Function to convert a list to np.float32 array.\n",
    "\n",
    "    Args:\n",
    "        long_list (list): A regular Python list with values that can be converted to floats.\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: A np.float32 array.\n",
    "\n",
    "    \"\"\"\n",
    "    np_array = (\n",
    "        np.zeros(\n",
    "            [len(max(long_list, key=lambda x: len(x))), len(long_list)],\n",
    "            dtype=np.float32,\n",
    "        )\n",
    "        - 1\n",
    "    )\n",
    "    for i, j in enumerate(long_list):\n",
    "        np_array[0 : len(j), i] = j\n",
    "\n",
    "    return np_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving output\n",
    "\n",
    "As MS hardware has continued to improve over the years, MS data has become more complex. To deal with this complexity, the MS community has already used many different [data formats](https://onlinelibrary.wiley.com/doi/full/10.1002/mas.21522) to store and access data. [HDF](https://www.hdfgroup.org/solutions/hdf5/) containers are one option, but they have not yet gained widespread support."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HDF containers in general\n",
    "\n",
    "In general, an HDF container can be viewed as a compressed folder with metadata (i.e., attributes) associated to each single subfolder or file (i.e., data arrays of various types and sizes) within this container. A container might, for instance have contents that look like, e.g.:\n",
    "```\n",
    "HDF_Container\n",
    "    {\n",
    "        meta_data_1: \"Some string\",\n",
    "        meta_data_2: 1234567890,\n",
    "        ...\n",
    "    }\n",
    "    array_1\n",
    "        {\n",
    "            meta_data_of_array1_1: \"Some other string\",\n",
    "            ...\n",
    "        },\n",
    "        100x2 int8\n",
    "    array_2\n",
    "        1000x2 float64\n",
    "    subfolder_1\n",
    "        {\n",
    "            meta_data_of_subfolder_1_1: \"Really any string of any length\",\n",
    "            ...\n",
    "        }\n",
    "        array_1_of_subfolder_1\n",
    "        subfolder_1_1\n",
    "        ...\n",
    "    subfolder_n\n",
    "    ...\n",
    "```\n",
    "\n",
    "\n",
    "A few of the advantages of HDF are, e.g.:\n",
    "\n",
    "* It has no upper limit with regards to file size\n",
    "* It can be used on disk without consuming RAM memory\n",
    "* It is fully portable on different machines and systems\n",
    "* It is capable of fast IO operations\n",
    "* It allows data to be structured very transparent while still providing a flexible way to store metadata\n",
    "* ...\n",
    "\n",
    "For these reasons, HDF containers have gained popularity in several scientific fields, including (astro)physics and geology. It is, therefore, no surprise that Python has excellent support for HDF containers. The two most used packages are `h5py` and `tables`, where the former has a generic API and the second is frequently used with `pandas` dataframes. An excellent viewer for HDF files is [HDF Compass](https://support.hdfgroup.org/projects/compass/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using HDF containers in AlphaPept\n",
    "\n",
    "AlphaPept uses the python `h5py` package to store MS data in HDF containers, inspired by the [ion_networks](https://github.com/swillems/ion_networks) repository.\n",
    "\n",
    "* First, a generic class is defined that will serve as an API for HDF containers. To ensure full transparency, this will include immutable metadata such as `creation time`, `original_file_name` and `version`.\n",
    "* The constructor of an HDF_File will be passed the `file_name` of an HDF container, an `is_read_only` flag, an `is_overwritable` flag and `is_new_file` flag.\n",
    "* To compare HDF_Files, several (magic) functions need to be defined.\n",
    "* Traceability and reproducibility are ensured by storing a `last_updated` and a `check` function to warn users about potential compatibility issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "import h5py\n",
    "import os\n",
    "import time\n",
    "from alphapept.__main__ import VERSION_NO\n",
    "\n",
    "\n",
    "class HDF_File(object):\n",
    "    '''\n",
    "    A generic class to store and retrieve on-disk\n",
    "    data with an HDF container.\n",
    "    '''\n",
    "\n",
    "    @property\n",
    "    def original_file_name(self):\n",
    "        return self.read(\n",
    "            attr_name=\"original_file_name\"\n",
    "        )  # See below for function definition\n",
    "\n",
    "    @property\n",
    "    def file_name(self):\n",
    "        return self.__file_name\n",
    "\n",
    "    @property\n",
    "    def directory(self):\n",
    "        return os.path.dirname(self.file_name)\n",
    "\n",
    "    @property\n",
    "    def creation_time(self):\n",
    "        return self.read(\n",
    "            attr_name=\"creation_time\"\n",
    "        )  # See below for function definition\n",
    "\n",
    "    @property\n",
    "    def last_updated(self):\n",
    "        return self.read(\n",
    "            attr_name=\"last_updated\"\n",
    "        )  # See below for function definition\n",
    "\n",
    "    @property\n",
    "    def version(self):\n",
    "        return self.read(\n",
    "            attr_name=\"version\"\n",
    "        )  # See below for function definition\n",
    "\n",
    "    @property\n",
    "    def is_read_only(self):\n",
    "        return self.__is_read_only\n",
    "\n",
    "    @property\n",
    "    def is_overwritable(self):\n",
    "        return self.__is_overwritable\n",
    "\n",
    "    def read(self):\n",
    "        pass\n",
    "\n",
    "    def write(self):\n",
    "        pass\n",
    "\n",
    "    def __init__(\n",
    "        self,\n",
    "        file_name: str,\n",
    "        is_read_only: bool = True,\n",
    "        is_new_file: bool = False,\n",
    "        is_overwritable: bool = False,\n",
    "    ):\n",
    "        \"\"\"Create/open a wrapper object to access HDF data.\n",
    "\n",
    "        Args:\n",
    "            file_name (str): The file_name of the HDF file.\n",
    "            is_read_only (bool): If True, the HDF file cannot be modified. Defaults to True.\n",
    "            is_new_file (bool): If True, an already existing file will be completely removed. Defaults to False.\n",
    "            is_overwritable (bool): If True, already existing arrays will be overwritten. If False, only new data can be appended. Defaults to False.\n",
    "\n",
    "        \"\"\"\n",
    "        self.__file_name = os.path.abspath(file_name)\n",
    "        if is_new_file:\n",
    "            is_read_only = False\n",
    "            if not os.path.exists(self.directory):\n",
    "                os.makedirs(self.directory)\n",
    "            with h5py.File(self.file_name, \"w\") as hdf_file:\n",
    "                current_time = time.asctime()\n",
    "                hdf_file.attrs[\"creation_time\"] = current_time\n",
    "                hdf_file.attrs[\"original_file_name\"] = self.__file_name\n",
    "                hdf_file.attrs[\"version\"] = VERSION_NO\n",
    "                hdf_file.attrs[\"last_updated\"] = current_time\n",
    "        else:\n",
    "            with h5py.File(self.file_name, \"r\") as hdf_file:\n",
    "                self.check()\n",
    "        if is_overwritable:\n",
    "            is_read_only = False\n",
    "        self.__is_read_only = is_read_only\n",
    "        self.__is_overwritable = is_overwritable\n",
    "\n",
    "    def __eq__(self, other):\n",
    "        return self.file_name == other.file_name\n",
    "\n",
    "    def __hash__(self):\n",
    "        return hash(self.file_name)\n",
    "\n",
    "    def __str__(self):\n",
    "        return f\"<HDF_File {self.file_name}>\"\n",
    "\n",
    "    def __repr__(self):\n",
    "        return str(self)\n",
    "\n",
    "    def check(\n",
    "        self,\n",
    "        version: bool = True,\n",
    "        file_name: bool = True,\n",
    "    ) -> list:\n",
    "        \"\"\"Check if the `version` or `file_name` of this HDF_File have changed.\n",
    "\n",
    "        Args:\n",
    "            version (bool): If False, do not check the version. Defaults to True.\n",
    "            file_name (bool): If False, do not check the file_name. Defaults to True.\n",
    "\n",
    "        Returns:\n",
    "            list: A list of warning messages stating any issues.\n",
    "\n",
    "        \"\"\"\n",
    "        warning_messages = []\n",
    "        if version:\n",
    "            current_version = VERSION_NO\n",
    "            creation_version = self.version\n",
    "            if creation_version != current_version:\n",
    "                warning_messages.append(\n",
    "                    f\"{self} was created with version \"\n",
    "                    f\"{creation_version} instead of {current_version}.\"\n",
    "                )\n",
    "        if file_name:\n",
    "            if self.file_name != self.original_file_name:\n",
    "                warning_messages.append(\n",
    "                    f\"The file name of {self} has been changed from\"\n",
    "                    f\"{self.original_file_name} to {self.file_name}.\"\n",
    "                )\n",
    "        return warning_messages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contents of HDF containers come in three variants:\n",
    "\n",
    "1. `Groups`: folders\n",
    "2. `Datasets`: arrays\n",
    "3. `Attributes`: metadata associated with individual datasets or groups (with the root folder also considered as a normal group)\n",
    "\n",
    "These contents can be accessed with `read` and `write` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "import pandas as pd\n",
    "from fastcore.foundation import patch\n",
    "\n",
    "\n",
    "@patch\n",
    "def read(\n",
    "    self: HDF_File,\n",
    "    group_name: str = None,\n",
    "    dataset_name: str = None,\n",
    "    attr_name: str = None,\n",
    "    return_dataset_shape: bool = False,\n",
    "    return_dataset_dtype: bool = False,\n",
    "    return_dataset_slice: slice = slice(None),\n",
    "    swmr: bool = False,\n",
    "):\n",
    "    \"\"\"Read contents of an HDF_File.\n",
    "\n",
    "    Args:\n",
    "        group_name (str): The group_name from where to read data.\n",
    "            If no `group_name` has been provided, read directly from the root group.\n",
    "            Defaults to None.\n",
    "        dataset_name (str): The dataset to read.\n",
    "            If no `dataset_name` has been provided, read directly from the group.\n",
    "            If the `dataset_name` refers to a group, it is assumed to be\n",
    "            pd.DataFrame and returned as such.\n",
    "            Defaults to None.\n",
    "        attr_name (str): The attribute to read.\n",
    "            If `attr_name` is not None, read the attribute value instead of the contents of a group or dataset.\n",
    "            If `attr_name` == \"\", read all attributes as a dict.\n",
    "            Defaults to None.\n",
    "        return_dataset_shape (bool): Do not read complete dataset to minimize RAM and IO usage.\n",
    "            Defaults to False.\n",
    "        return_dataset_dtype (bool): Do not read complete dataset to minimize RAM and IO usage.\n",
    "            Defaults to False.\n",
    "        return_dataset_slice (slice): Do not read complete dataset to minimize RAM and IO usage.\n",
    "            Defaults to slice(None).\n",
    "        swmr (bool): Use swmr mode to read data. Defaults to False.\n",
    "\n",
    "    Returns:\n",
    "        type: Depending on what is requested, a dict, value, np.ndarray or pd.dataframe is returned.\n",
    "\n",
    "    Raises:\n",
    "        KeyError: When the group_name does not exist.\n",
    "        KeyError: When the attr_name does not exist in the group or dataset.\n",
    "        KeyError: When the dataset_name does not exist in the group.\n",
    "        ValueError: When the requested dataset is not a np.ndarray or pd.dataframe.\n",
    "\n",
    "    \"\"\"\n",
    "    with h5py.File(self.file_name, \"r\", swmr=swmr) as hdf_file:\n",
    "        if group_name is None:\n",
    "            group = hdf_file\n",
    "            group_name = \"/\"\n",
    "        else:\n",
    "            try:\n",
    "                group = hdf_file[group_name]\n",
    "            except KeyError as k:\n",
    "                raise KeyError(\n",
    "                    f\"Group {group_name} does not exist in {self}. Error {k}\"\n",
    "                )\n",
    "        if dataset_name is None:\n",
    "            if attr_name is None:\n",
    "                return sorted(group)\n",
    "            elif attr_name != \"\":\n",
    "                try:\n",
    "                    return group.attrs[attr_name]\n",
    "                except KeyError:\n",
    "                    raise KeyError(\n",
    "                        f\"Attribute {attr_name} does not exist for \"\n",
    "                        f\"group {group_name} of {self}.\"\n",
    "                    )\n",
    "            else:\n",
    "                return dict(group.attrs)\n",
    "        else:\n",
    "            try:\n",
    "                dataset = group[dataset_name]\n",
    "            except KeyError:\n",
    "                raise KeyError(\n",
    "                    f\"Dataset {dataset_name} does not exist for \"\n",
    "                    f\"group {group_name} of {self}.\"\n",
    "                )\n",
    "            if attr_name is None:\n",
    "                if isinstance(dataset, h5py.Dataset):\n",
    "                    if return_dataset_shape:\n",
    "                        return dataset.shape\n",
    "                    elif return_dataset_dtype:\n",
    "                        return dataset.dtype\n",
    "                    else:\n",
    "                        array = dataset[return_dataset_slice]\n",
    "                        # TODO: This assumes any object array is a string array\n",
    "                        if array.dtype == object:\n",
    "                            array = array.astype(str)\n",
    "                        return array\n",
    "                elif dataset.attrs[\"is_pd_dataframe\"]:\n",
    "                    if return_dataset_shape:\n",
    "                        columns = list(dataset)\n",
    "                        return (\n",
    "                            len(dataset[columns[0]]),\n",
    "                            len(columns)\n",
    "                        )\n",
    "                    elif return_dataset_dtype:\n",
    "                        return [\n",
    "                            dataset[column].dtype for column in sorted(\n",
    "                                dataset\n",
    "                            )\n",
    "                        ]\n",
    "                    else:\n",
    "                        df = pd.DataFrame(\n",
    "                            {\n",
    "                                column: dataset[column][\n",
    "                                    return_dataset_slice\n",
    "                                ] for column in sorted(dataset)\n",
    "                            }\n",
    "                        )\n",
    "                        # TODO: This assumes any object array is a string array\n",
    "                        for column in dataset:\n",
    "                            if df[column].dtype == object:\n",
    "                                df[column] = df[column].apply(\n",
    "                                    lambda x: x if isinstance(x, str) else x.decode('UTF-8')\n",
    "                                )\n",
    "                        return df\n",
    "                else:\n",
    "                    raise ValueError(\n",
    "                        f\"{dataset_name} is not a valid dataset in \"\n",
    "                        f\"group {group_name} of {self}.\"\n",
    "                    )\n",
    "            elif attr_name != \"\":\n",
    "                try:\n",
    "                    return dataset.attrs[attr_name]\n",
    "                except KeyError:\n",
    "                    raise KeyError(\n",
    "                        f\"Attribute {attr_name} does not exist for \"\n",
    "                        f\"dataset {dataset_name} of group \"\n",
    "                        f\"{group_name} of {self}.\"\n",
    "                    )\n",
    "            else:\n",
    "                return dict(dataset.attrs)\n",
    "\n",
    "\n",
    "@patch\n",
    "def write(\n",
    "    self:HDF_File,\n",
    "    value,\n",
    "    group_name:str=None,\n",
    "    dataset_name:str=None,\n",
    "    attr_name:str=None,\n",
    "    overwrite:bool=None,\n",
    "    dataset_compression:str=None,\n",
    "    swmr:bool=False,\n",
    ") -> None:\n",
    "    \"\"\"Write a `value` to an HDF_File.\n",
    "\n",
    "    Args:\n",
    "        value (type): The name of the data to write.\n",
    "            If the `value` is pd.DataFrame, a `dataset_name` must be provided.\n",
    "        group_name (str): The group where to write data.\n",
    "            If no `group_name` is provided, write directly to the root group.\n",
    "            Defaults to None.\n",
    "        dataset_name (str): If no `dataset_name` is provided,\n",
    "            create a new group with `value` as name.\n",
    "            The dataset where to write data. Defaults to None.\n",
    "        attr_name (str): The attr where to write data. Defaults to None.\n",
    "        overwrite (bool): Overwrite pre-existing data and truncate existing groups.\n",
    "            If the False, ignore the is_overwritable flag of this HDF_File.\n",
    "            Defaults to None.\n",
    "        dataset_compression (str): The compression type to use for datasets.\n",
    "            Defaults to None.\n",
    "        swmr (bool): Open files in swmr mode. Defaults to False.\n",
    "\n",
    "    Raises:\n",
    "        IOError: When the object is read-only.\n",
    "        KeyError: When the group_name or attr_name does not exist.\n",
    "        ValueError: When trying to overwrite something while overwiting is disabled.\n",
    "\n",
    "    \"\"\"\n",
    "    if self.is_read_only:\n",
    "        raise IOError(\n",
    "            f\"Trying to write to {self}, which is read_only.\"\n",
    "        )\n",
    "    if overwrite is None:\n",
    "        overwrite = self.is_overwritable\n",
    "    with h5py.File(self.file_name, \"a\", swmr=swmr) as hdf_file:\n",
    "\n",
    "        if group_name is None:\n",
    "            group = hdf_file\n",
    "            group_name = \"/\"\n",
    "        else:\n",
    "            try:\n",
    "                group = hdf_file[group_name]\n",
    "            except KeyError:\n",
    "                raise KeyError(\n",
    "                    f\"Group {group_name} does not exist in {self}.\"\n",
    "                )\n",
    "        if dataset_name is None:\n",
    "            if attr_name is None:\n",
    "                if value in group:\n",
    "                    if overwrite:\n",
    "                        del group[value]\n",
    "                    else:\n",
    "                        raise ValueError(\n",
    "                            f\"New group {value} already exists in group \"\n",
    "                            f\"{group_name} of {self}.\"\n",
    "                        )\n",
    "                group.create_group(value)\n",
    "            else:\n",
    "                if (attr_name in group.attrs) and not overwrite:\n",
    "                    raise ValueError(\n",
    "                        f\"Attribute {attr_name} already exists in group \"\n",
    "                        f\"{group_name} of {self}.\"\n",
    "                    )\n",
    "                try:\n",
    "                    group.attrs[attr_name] = value\n",
    "                except TypeError:\n",
    "                    group.attrs[attr_name] = str(value)\n",
    "        else:\n",
    "            if attr_name is None:\n",
    "                if dataset_name in group:\n",
    "                    if overwrite:\n",
    "                        del group[dataset_name]\n",
    "                    else:\n",
    "                        raise ValueError(\n",
    "                            f\"Dataset {dataset_name} already exists in group \"\n",
    "                            f\"{group_name} of {self}.\"\n",
    "                        )\n",
    "                if isinstance(value, pd.core.frame.DataFrame):\n",
    "                    new_group_name = f\"{group_name}/{dataset_name}\"\n",
    "                    self.write(\n",
    "                        dataset_name,\n",
    "                        group_name=group_name,\n",
    "                        overwrite=overwrite,\n",
    "                    )\n",
    "                    self.write(\n",
    "                        True,\n",
    "                        group_name=new_group_name,\n",
    "                        attr_name=\"is_pd_dataframe\",\n",
    "                        overwrite=overwrite,\n",
    "                    )\n",
    "                    for column in value.columns:\n",
    "                        self.write(\n",
    "                            value[column].values,\n",
    "                            group_name=new_group_name,\n",
    "                            dataset_name=column,\n",
    "                            overwrite=overwrite,\n",
    "                            dataset_compression=dataset_compression,\n",
    "                        )\n",
    "                else:\n",
    "                    dtype = value.dtype\n",
    "                    if value.dtype == np.dtype('O'):\n",
    "                        dtype = h5py.string_dtype()\n",
    "                    try:\n",
    "                        hdf_dataset = group.create_dataset(\n",
    "                            dataset_name,\n",
    "                            data=value,\n",
    "                            compression=dataset_compression,\n",
    "                            dtype=dtype\n",
    "                        )\n",
    "                    except TypeError:\n",
    "                        # TODO\n",
    "                        # print(f\"Cannot save array {value} to HDF, skipping it...\")\n",
    "                        pass\n",
    "            else:\n",
    "                try:\n",
    "                    dataset = group[dataset_name]\n",
    "                except KeyError:\n",
    "                    raise KeyError(\n",
    "                        f\"Dataset {dataset_name} does not exist for \"\n",
    "                        f\"group {group_name} of {self}.\"\n",
    "                    )\n",
    "                if (attr_name in dataset.attrs) and not overwrite:\n",
    "                    raise ValueError(\n",
    "                        f\"Attribute {attr_name} already exists in \"\n",
    "                        f\"dataset {dataset_name} of group \"\n",
    "                        f\"{group_name} of {self}.\"\n",
    "                    )\n",
    "                try:\n",
    "                    dataset.attrs[attr_name] = value\n",
    "                except TypeError:\n",
    "                    dataset.attrs[attr_name] = str(value) # e.g. dicts\n",
    "        hdf_file.attrs[\"last_updated\"] = time.asctime()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unit tests for this generic HDF class include:\n",
    "\n",
    "* Creation and truncation of files with various access.\n",
    "* Writing and reading data from the container."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "\n",
    "import numpy as np\n",
    "import unittest\n",
    "import os\n",
    "import alphapept.io\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "def define_new_test_files(test_folder):\n",
    "    test_file_names = [\n",
    "        os.path.abspath(\n",
    "            os.path.join(test_folder, f\"{file_name}.hdf\")\n",
    "        ) for file_name in [\n",
    "            \"test0\",\n",
    "            \"test1\",\n",
    "            \"test2\",\n",
    "        ]\n",
    "    ]\n",
    "    for file_name in test_file_names:\n",
    "        if os.path.isfile(file_name):\n",
    "            os.remove(file_name)\n",
    "    return test_file_names\n",
    "\n",
    "def test_hdf_file_creation(test_folder):\n",
    "    test_file_names = define_new_test_files(test_folder)\n",
    "    try:\n",
    "        f0 = HDF_File(test_file_names[0])\n",
    "    except OSError:\n",
    "        assert True\n",
    "    else:\n",
    "        assert False, \"Non-existing file should raise an error\"\n",
    "    f0 = HDF_File(test_file_names[0], is_new_file=True)\n",
    "    assert f0.is_read_only == False, \"New files should never be read-only\"\n",
    "    del f0\n",
    "    try:\n",
    "        f0 = HDF_File(test_file_names[0])\n",
    "    except OSError:\n",
    "        assert False, \"Newly created file should exist on disk\"\n",
    "    else:\n",
    "        assert True\n",
    "    assert f0.is_read_only == True, \"Existing files should be read-only\"\n",
    "    assert f0.file_name == test_file_names[0], \"File name should match given file name\"\n",
    "    assert f0.original_file_name == test_file_names[0], \"Original file name should match given file name\"\n",
    "    assert f0.version == VERSION_NO, \"Versions should match\"\n",
    "    assert str(f0) == f\"<HDF_File {test_file_names[0]}>\", \"File name should match\"\n",
    "    del f0\n",
    "    os.rename(test_file_names[0], test_file_names[1])\n",
    "    f1 = HDF_File(test_file_names[1])\n",
    "    assert f1.original_file_name != test_file_names[1], \"Original file name should not match given file name\"\n",
    "    assert f1.directory == os.path.abspath(test_folder), \"Directory should match\"\n",
    "    f0 = HDF_File(test_file_names[0], is_new_file=True)\n",
    "    assert f0 != f1, \"Different file names should be different HDF_Files\"\n",
    "    del f0\n",
    "    f1_copy = HDF_File(test_file_names[1], is_read_only=False)\n",
    "    assert f1 == f1_copy, \"Same file names should be same HDF_Files\"\n",
    "    assert f1_copy.is_read_only == False, \"File should not be read-only\"\n",
    "\n",
    "\n",
    "def test_hdf_file_read_and_write(test_folder):\n",
    "    test_file_names = define_new_test_files(test_folder)\n",
    "    f0 = HDF_File(test_file_names[0], is_new_file=True)\n",
    "    try:\n",
    "        f0.read(group_name=\"subgroup\")\n",
    "    except KeyError:\n",
    "        assert True\n",
    "    else:\n",
    "        assert False, \"subgroup should not exist\"\n",
    "    f0.write(\"subgroup\")\n",
    "    try:\n",
    "        f0.read(group_name=\"subgroup\")\n",
    "    except KeyError:\n",
    "        assert False, \"Subgroup should exist\"\n",
    "    else:\n",
    "        assert True\n",
    "    z = np.random.random((100, 4))\n",
    "    f0.write(z, group_name=\"subgroup\", dataset_name=\"random\")\n",
    "    f0_copy = HDF_File(test_file_names[0])\n",
    "    assert f0_copy.read(\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        return_dataset_shape=True\n",
    "    ) == z.shape, \"Shape of dataset is not correct\"\n",
    "    assert f0_copy.read(\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        return_dataset_dtype=True\n",
    "    ) == z.dtype, \"Type of dataset is not correct\"\n",
    "    assert np.all(\n",
    "        f0_copy.read(dataset_name=\"random\", group_name=\"subgroup\") == z\n",
    "    ), \"Contents of dataset are not correct\"\n",
    "    try:\n",
    "        f0.write(z, group_name=\"subgroup\", dataset_name=\"random\")\n",
    "    except ValueError:\n",
    "        assert True\n",
    "    else:\n",
    "        assert False, \"Should not overwrite dataset\"\n",
    "    try:\n",
    "        f0.write(z, group_name=\"subgroup\", dataset_name=\"random\", overwrite=True)\n",
    "    except ValueError:\n",
    "        assert False, \"Should be able to overwrite dataset\"\n",
    "    else:\n",
    "        assert True\n",
    "    f0.write(\n",
    "        4,\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"numeric_attr\",\n",
    "    )\n",
    "    assert f0.read(\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"numeric_attr\",\n",
    "    ) == 4, \"Attr shoud match\"\n",
    "    f0.write(\n",
    "        \"test\",\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"string_attr\",\n",
    "    )\n",
    "    assert f0.read(\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"string_attr\",\n",
    "    ) == \"test\", \"String attr shoud match\"\n",
    "    f0.write(\n",
    "        list(range(5)),\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"int_list_attr\",\n",
    "    )\n",
    "    assert np.all(\n",
    "        f0.read(\n",
    "            dataset_name=\"random\",\n",
    "            group_name=\"subgroup\",\n",
    "            attr_name=\"int_list_attr\",\n",
    "        ) == list(range(5))\n",
    "    ), \"Attr shoud match\"\n",
    "    mixed_list =  [\"test\", \"mixed\", 2, 4.9]\n",
    "    try:\n",
    "        f0.write(\n",
    "            mixed_list,\n",
    "            dataset_name=\"random\",\n",
    "            group_name=\"subgroup\",\n",
    "            attr_name=\"int_list_attr\",\n",
    "        )\n",
    "    except ValueError:\n",
    "        assert True\n",
    "    else:\n",
    "        assert False, \"Should not be able to overwrite attr\"\n",
    "    f0.write(\n",
    "        mixed_list,\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"int_list_attr\",\n",
    "        overwrite=True\n",
    "    )\n",
    "#     assert np.all(\n",
    "#         f0.read(\n",
    "#             dataset_name=\"random\",\n",
    "#             group_name=\"subgroup\",\n",
    "#             attr_name=\"int_list_attr\",\n",
    "#         ) == mixed_list\n",
    "#     ), \"Attr of mixed list should match\"\n",
    "    f0.write(\n",
    "        {\"t\": 1},\n",
    "        dataset_name=\"random\",\n",
    "        group_name=\"subgroup\",\n",
    "        attr_name=\"dict_attr\",\n",
    "        overwrite=True\n",
    "    )\n",
    "    assert np.all(\n",
    "        f0.read(\n",
    "            dataset_name=\"random\",\n",
    "            group_name=\"subgroup\",\n",
    "            attr_name=\"dict_attr\",\n",
    "        ) == str({\"t\": 1})\n",
    "    ), \"Attr shoud match\"\n",
    "    \n",
    "\n",
    "def test_hdf_file_data_frames(test_folder):\n",
    "    test_file_names = define_new_test_files(test_folder)\n",
    "    f0 = HDF_File(test_file_names[0], is_new_file=True)\n",
    "    df = pd.DataFrame(\n",
    "        {\n",
    "            \"col1\": np.arange(10, dtype=np.float16) / 2,\n",
    "            \"col2\": np.arange(10),\n",
    "        }\n",
    "    )\n",
    "    f0.write(df, dataset_name=\"df\")\n",
    "    z = f0.read(dataset_name=\"df\")\n",
    "    assert z.equals(df)\n",
    "    \n",
    "test_hdf_file_creation(test_folder=\"tmp\")\n",
    "test_hdf_file_read_and_write(test_folder=\"tmp\")\n",
    "test_hdf_file_data_frames(test_folder=\"tmp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The final `.ms_data.hdf` file\n",
    "\n",
    "Based on the generic HDF_File, a subclass that acts as an MS data container can be implemented. This class should contain all (centroided) fragment ions, all their coordinates, and all the metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "class MS_Data_File(HDF_File):\n",
    "    \"\"\" A class to store and retrieve on-disk MS data with an HDF container.\"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single generic function should allow to read raw data and store spectra. Different arguments allow different vendor formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "@patch\n",
    "def import_raw_DDA_data(\n",
    "    self:MS_Data_File,\n",
    "    file_name:str,\n",
    "    n_most_abundant:int=-1,\n",
    "    callback:callable=None,\n",
    "    query_data:dict=None,\n",
    "    vendor:str=None\n",
    ") -> None:\n",
    "    \"\"\"Load centroided data and save it to this object.\n",
    "\n",
    "    Args:\n",
    "        file_name (str): The file name with raw data (Thermo, Bruker or mzml).\n",
    "        n_most_abundant (int): The maximum number of peaks to retain per MS2 spectrum.\\\n",
    "            Defaults to -1.\n",
    "        callback (callable): A function that accepts a float between 0 and 1 as progress. Defaults to None.\n",
    "        query_data (dict):\n",
    "            A dictionary with raw data.\n",
    "            If this is not None, data will only be saved and not imported.\n",
    "            Defaults to None.\n",
    "        vendor (str): The vendor name, must be Thermo or Bruker if provided.\n",
    "            Defaults to None.\n",
    "\n",
    "    \"\"\"\n",
    "    base, ext = os.path.splitext(file_name)\n",
    "    if query_data is None:\n",
    "        query_data, vendor, acquisition_date_time = _read_DDA_query_data(\n",
    "            file_name,\n",
    "            n_most_abundant=n_most_abundant,\n",
    "            callback=callback\n",
    "        )\n",
    "    self._save_DDA_query_data(query_data, vendor, acquisition_date_time)\n",
    "    \n",
    "    \n",
    "def index_ragged_list(ragged_list: list)  -> np.ndarray:\n",
    "    \"\"\"Create lookup indices for a list of arrays for concatenation.\n",
    "\n",
    "    Args:\n",
    "        value (list): Input list of arrays.\n",
    "\n",
    "    Returns:\n",
    "        indices: A numpy array with indices.\n",
    "    \"\"\"\n",
    "    indices = np.zeros(len(ragged_list) + 1, np.int64)\n",
    "    indices[1:] = [len(i) for i in ragged_list]\n",
    "    indices = np.cumsum(indices)\n",
    "    \n",
    "    return indices\n",
    "\n",
    "def _read_DDA_query_data(\n",
    "    file_name:str,\n",
    "    n_most_abundant:int=-1,\n",
    "    callback:callable=None\n",
    ") -> tuple:\n",
    "    \"\"\"Read raw data and return as query dictionary.\n",
    "\n",
    "    Args:\n",
    "        file_name (str): The file name with raw data (Thermo, Bruker or mzml).\n",
    "        n_most_abundant (int): The maximum number of peaks to retain per MS2 spectrum.\\\n",
    "            Defaults to -1.\n",
    "        callback (callable): A function that accepts a float between 0 and 1 as progress. Defaults to None.\n",
    "\n",
    "    Returns:\n",
    "        tuple: A tuple with (query_data, vendor, acquisition_date_time).\n",
    "\n",
    "    Raises:\n",
    "        NotImplementedError: If the raw data is no Bruker, Thermo or mzml.\n",
    "\n",
    "    \"\"\"\n",
    "    base, ext = os.path.splitext(file_name)\n",
    "    if ext.lower() == '.raw':\n",
    "        if os.path.isdir(file_name):\n",
    "            vendor = \"Waters\"\n",
    "            acquisition_date_time = None\n",
    "            raise NotImplementedError(\n",
    "                f'File extension {ext} indicates Waters, which is not implemented.'\n",
    "            )\n",
    "        else:\n",
    "            vendor = \"Thermo\"\n",
    "            logging.info(f'File {base} has extension {ext} - converting from {vendor}.')\n",
    "            query_data, acquisition_date_time = load_thermo_raw(\n",
    "                file_name,\n",
    "                n_most_abundant,\n",
    "                callback=callback,\n",
    "            )\n",
    "    elif ext.lower() == '.wiff':\n",
    "        vendor = \"SciEx\"\n",
    "        logging.info(f'File {base} has extension {ext} - converting from {vendor}.')\n",
    "        query_data, acquisition_date_time = load_sciex_raw(\n",
    "            file_name,\n",
    "            n_most_abundant,\n",
    "            callback=callback,\n",
    "        )\n",
    "    elif ext.lower() == '.d':\n",
    "        vendor = \"Bruker\"\n",
    "        logging.info(f'File {base} has extension {ext} - converting from {vendor}.')\n",
    "        query_data, acquisition_date_time = load_bruker_raw(\n",
    "            file_name,\n",
    "            n_most_abundant,\n",
    "            callback=callback,\n",
    "        )\n",
    "    elif ext.lower() == '.mzml':\n",
    "        logging.info(f'File {base} has extension {ext} - converting from {ext[1:]}.')\n",
    "        query_data, acquisition_date_time, vendor = load_mzml_data(\n",
    "            file_name,\n",
    "            n_most_abundant,\n",
    "            callback=callback,\n",
    "        )\n",
    "    else:\n",
    "        raise NotImplementedError(f'File extension {ext} not understood.')\n",
    "    logging.info(\n",
    "        f'File conversion complete. Extracted {len(query_data[\"prec_mass_list2\"]):,} precursors.'\n",
    "    )\n",
    "    return query_data, vendor, acquisition_date_time\n",
    "\n",
    "@patch\n",
    "def _save_DDA_query_data(\n",
    "    self:MS_Data_File,\n",
    "    query_data:dict,\n",
    "    vendor:str,\n",
    "    acquisition_date_time:str,\n",
    "    overwrite=False\n",
    "):\n",
    "#     if vendor == \"Bruker\":\n",
    "#         raise NotImplementedError(\"Unclear what are ms1 and ms2 attributes for bruker\")\n",
    "    if \"Raw\" not in self.read():\n",
    "        self.write(\"Raw\")\n",
    "    self.write(vendor, group_name=\"Raw\", attr_name=\"vendor\")\n",
    "    self.write(acquisition_date_time, group_name=\"Raw\", attr_name=\"acquisition_date_time\")\n",
    "    if \"MS1_scans\" not in self.read(group_name=\"Raw\"):\n",
    "        self.write(\"MS1_scans\", group_name=\"Raw\")\n",
    "    if \"MS2_scans\" not in self.read(group_name=\"Raw\"):\n",
    "        self.write(\"MS2_scans\", group_name=\"Raw\")\n",
    "    for key, value in query_data.items():\n",
    "        if key.endswith(\"1\"):\n",
    "#             TODO: Weak check for ms2, imporve to _ms1 if consistency in naming is guaranteed\n",
    "            if key == \"mass_list_ms1\":\n",
    "                indices = np.zeros(len(value) + 1, np.int64)\n",
    "                indices[1:] = [len(i) for i in value]\n",
    "                indices = np.cumsum(indices)\n",
    "                self.write(\n",
    "                    indices,\n",
    "                    dataset_name=\"indices_ms1\",\n",
    "                    group_name=f\"Raw/MS1_scans\"\n",
    "                )\n",
    "                value = np.concatenate(value)\n",
    "            elif key == \"int_list_ms1\":\n",
    "                value = np.concatenate(value)\n",
    "            self.write(\n",
    "                value,\n",
    "#                 TODO: key should be trimmed: xxx_ms1 should just be e.g. xxx\n",
    "                dataset_name=key,\n",
    "                group_name=f\"Raw/MS1_scans\"\n",
    "            )\n",
    "        elif key.endswith(\"2\"):\n",
    "#             TODO: Weak check for ms2, imporve to _ms2 if consistency in naming is guaranteed\n",
    "            if key == \"mass_list_ms2\":\n",
    "                indices = np.zeros(len(value) + 1, np.int64)\n",
    "                indices[1:] = [len(i) for i in value]\n",
    "                indices = np.cumsum(indices)\n",
    "                self.write(\n",
    "                    indices,\n",
    "                    dataset_name=\"indices_ms2\",\n",
    "                    group_name=f\"Raw/MS2_scans\"\n",
    "                )\n",
    "                value = np.concatenate(value)\n",
    "            elif key == \"int_list_ms2\":\n",
    "                value = np.concatenate(value)\n",
    "            elif key == \"alphatims_spectrum_indptr_ms2\":\n",
    "                key = \"indices_ms2\"\n",
    "            elif key == \"alphatims_spectrum_mz_values_ms2\":\n",
    "                key = \"mass_list_ms2\"\n",
    "            elif key == \"alphatims_spectrum_intensity_values_ms2\":\n",
    "                key = \"int_list_ms2\"\n",
    "            self.write(\n",
    "                value,\n",
    "#                 TODO: key should be trimmed: xxx_ms2 should just be e.g. xxx\n",
    "                dataset_name=key,\n",
    "                group_name=f\"Raw/MS2_scans\"\n",
    "            )\n",
    "        else:\n",
    "            raise KeyError(\"Unspecified scan type\")\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testing of the MS_Data_File container includes reading and writing from different file formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "\n",
    "def test_ms_data_file_import_thermo(file_name):\n",
    "    local_file_name = os.path.basename(file_name)\n",
    "    output_path = os.path.dirname(file_name)\n",
    "    base_file_name, ext = os.path.splitext(local_file_name)\n",
    "    output_file_name = os.path.join(output_path, base_file_name+\".ms_data.hdf\")\n",
    "    ms_data_file = MS_Data_File(\n",
    "        output_file_name,\n",
    "        is_new_file=True\n",
    "    )\n",
    "    ms_data_file.import_raw_DDA_data(\n",
    "        file_name,\n",
    "        n_most_abundant = settings[\"raw\"][\"n_most_abundant\"]\n",
    "    )\n",
    "\n",
    "\n",
    "# print(time.asctime())\n",
    "# test_ms_data_file_import_thermo(\n",
    "#     \"/Users/swillems/Documents/sandbox/alphapept_projects/09-07-18_EcoliSpikeIn_1xF1R1.raw\"\n",
    "# )\n",
    "# print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While that HDF data structure could be used directly, it is often easier to read it and return a `query_data` dictionary similar to those that are returned by the readers of `Thermo`, `Bruker`, `mzML` and `mzXML` raw data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "@patch\n",
    "def read_DDA_query_data(\n",
    "    self:MS_Data_File,\n",
    "    calibrated_fragments:bool=False,\n",
    "    force_recalibrate:bool=False,\n",
    "    swmr:bool=False,\n",
    "    **kwargs\n",
    ") -> dict:\n",
    "    \"\"\"Read query data from this ms_data object and return it as a query_dict.\n",
    "\n",
    "    Args:\n",
    "        calibrated_fragments (bool): If True, calibrated fragments are retrieved.\n",
    "            Calibration offsets can already be present in the ms_data or recalculated.\n",
    "            Defaults to False.\n",
    "        force_recalibrate (bool): If calibrated fragments is True,\n",
    "            recalibrate mzs values even if a recalibration is already provided.\n",
    "            Defaults to False.\n",
    "        swmr (bool): Open the file in swmr mode. Defaults to False.\n",
    "        **kwargs (type): Can contain a database file name that was used for recalibration.\n",
    "\n",
    "    Returns:\n",
    "        dict: A query_dict with data for MS1 and MS2 scans.\n",
    "\n",
    "    \"\"\"\n",
    "    query_data = {}\n",
    "    for dataset_name in self.read(group_name=\"Raw/MS1_scans\"):\n",
    "        values = self.read(\n",
    "            dataset_name=dataset_name,\n",
    "            group_name=\"Raw/MS1_scans\",\n",
    "            swmr=swmr,\n",
    "        )\n",
    "        query_data[dataset_name] = values\n",
    "    for dataset_name in self.read(group_name=\"Raw/MS2_scans\"):\n",
    "        values = self.read(\n",
    "            dataset_name=dataset_name,\n",
    "            group_name=\"Raw/MS2_scans\",\n",
    "            swmr=swmr\n",
    "        )\n",
    "        query_data[dataset_name] = values\n",
    "#     indices_ms1 = query_data[\"indices_ms1\"]\n",
    "#     mz_ms1 = query_data[\"mass_list_ms1\"]\n",
    "#     query_data[\"mass_list_ms1\"] = np.array(\n",
    "#         [mz_ms1[s:e] for s,e in zip(indices_ms1[:-1], indices_ms1[1:])]\n",
    "#     )\n",
    "#     int_ms1 = query_data[\"int_list_ms1\"]\n",
    "#     query_data[\"int_list_ms1\"] = np.array(\n",
    "#         [int_ms1[s:e] for s,e in zip(indices_ms1[:-1], indices_ms1[1:])]\n",
    "#     )\n",
    "    indices_ms2 = query_data[\"indices_ms2\"]\n",
    "#     mz_ms2 = query_data[\"mass_list_ms2\"]\n",
    "#     query_data[\"mass_list_ms2\"] = np.array(\n",
    "#         [mz_ms2[s:e] for s,e in zip(indices_ms2[:-1], indices_ms2[1:])]\n",
    "#     )\n",
    "#     int_ms2 = query_data[\"int_list_ms2\"]\n",
    "#     query_data[\"int_list_ms2\"] = np.array(\n",
    "#         [int_ms2[s:e] for s,e in zip(indices_ms2[:-1], indices_ms2[1:])]\n",
    "#     )\n",
    "    if self.read(attr_name=\"vendor\", group_name=\"Raw\") == \"Bruker\":\n",
    "        query_data[\"mobility\"] = query_data[\"mobility2\"]\n",
    "        query_data[\"prec_id\"] = query_data[\"prec_id2\"]\n",
    "    if calibrated_fragments:\n",
    "        if (\"corrected_fragment_mzs\" not in self.read()) or force_recalibrate:\n",
    "#         if True:\n",
    "            logging.info(\"Calibrating fragments\")\n",
    "            import alphapept.recalibration\n",
    "            alphapept.recalibration.calibrate_fragments(\n",
    "                kwargs[\"database_file_name\"],\n",
    "                self.file_name,\n",
    "            )\n",
    "        query_data[\"mass_list_ms2\"] *= (\n",
    "            1 - self.read(\n",
    "                dataset_name=\"corrected_fragment_mzs\", swmr=swmr\n",
    "            ) / 10**6\n",
    "        )\n",
    "    return query_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "\n",
    "def test_get_query_datafrom_thermo_ms_file(file_name):\n",
    "    local_file_name = os.path.basename(file_name)\n",
    "    output_path = os.path.dirname(file_name)\n",
    "    base_file_name, ext = os.path.splitext(local_file_name)\n",
    "    output_file_name = os.path.join(output_path, base_file_name+\".ms_data.hdf\")\n",
    "    ms_data_file = MS_Data_File(\n",
    "        output_file_name,\n",
    "    )\n",
    "    return ms_data_file.read_DDA_query_data()\n",
    "\n",
    "\n",
    "# print(time.asctime())\n",
    "# qd = test_get_query_datafrom_thermo_ms_file(\n",
    "#     \"/Users/swillems/Documents/sandbox/alphapept_projects/09-07-18_EcoliSpikeIn_1xF1R1.raw\"\n",
    "# )\n",
    "# print(time.asctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Usage in workflows\n",
    "\n",
    "To use all the above functionality from a workflow with several parameters, the following functions are defined. These functions also allow parallel processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "\n",
    "def raw_conversion(\n",
    "    to_process: dict,\n",
    "    callback: callable = None,\n",
    "    parallel:bool = False\n",
    ") -> bool:\n",
    "    \"\"\"Wrapper function to convert raw to ms_data_file hdf.\n",
    "\n",
    "    Args:\n",
    "        to_process (dict): A dictionary with settings indicating which files are to be processed and how.\n",
    "        callback (callable): A function that accepts a float between 0 and 1 as progress. Defaults to None.\n",
    "        parallel (bool): If True, process multiple files in parallel.\n",
    "            This is not implemented yet!\n",
    "            Defaults to False.\n",
    "\n",
    "    Returns:\n",
    "        bool: True if and only if the conversion was succesful.\n",
    "\n",
    "    \"\"\"\n",
    "    index, settings = to_process\n",
    "    file_name = settings['experiment']['file_paths'][index]\n",
    "    try:\n",
    "        local_file_name = os.path.basename(file_name)\n",
    "        output_path = os.path.dirname(file_name)\n",
    "        base_file_name, ext = os.path.splitext(local_file_name)\n",
    "        output_file_name = os.path.join(output_path, base_file_name+\".ms_data.hdf\")\n",
    "\n",
    "        if not os.path.isfile(output_file_name):\n",
    "            ms_data_file = MS_Data_File(\n",
    "                output_file_name,\n",
    "                is_new_file=True\n",
    "            )\n",
    "            ms_data_file.import_raw_DDA_data(\n",
    "                file_name,\n",
    "                n_most_abundant = settings[\"raw\"][\"n_most_abundant\"]\n",
    "            )\n",
    "\n",
    "        logging.info(f'File conversion of file {file_name} complete.')\n",
    "        return True\n",
    "    except Exception as e:\n",
    "        logging.error(f'File conversion of file {file_name} failed. Exception {e}')\n",
    "        return f\"{e}\" #Can't return exception object, cast as string\n",
    "    return True\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following commands are bookkeeping to make sure this and other notebooks are properly parsed to python modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#|hide\n",
    "import nbdev; nbdev.nbdev_export()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}