project/10_0_ald_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "88e88633-9d3f-4751-a0aa-a508fb015203",
   "metadata": {
    "tags": []
   },
   "source": [
    "# ALD Study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "837e3205-4974-4bd4-bb39-288408fcdddc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import yaml\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pimmslearn\n",
    "\n",
    "logger = pimmslearn.logging.setup_nb_logger()\n",
    "\n",
    "pd.options.display.max_columns = 50\n",
    "pd.options.display.max_rows = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a1961d4-16d0-4eec-81ad-cc609847e0a9",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "folder_data = Path('data/ALD_study/')\n",
    "folder_data_out = folder_data / 'processed'\n",
    "folder_data_out.mkdir(parents=True, exist_ok=True)\n",
    "folder_run = Path('runs/appl_ald_data')\n",
    "folder_run.mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "print(*(folder_data.iterdir()), sep='\\n')\n",
    "\n",
    "fnames = dict(\n",
    "    plasma_proteinGroups=folder_data / 'Protein_ALDupgrade_Report.csv',\n",
    "    plasma_aggPeptides=folder_data / 'ald_proteome_spectronaut.tsv',\n",
    "    liver_proteinGroups=folder_data / 'Protein_20200221_121354_20200218_ALD_LiverTissue_PlateS1_Atlaslib_Report.csv',\n",
    "    liver_aggPeptides=folder_data / 'Peptide_20220819_100847_20200218_ALD_LiverTissue_PlateS1_Atlaslib_Report.csv',\n",
    "    annotations=folder_data / 'ald_experiment_annotations.csv',\n",
    "    clinic=folder_data / 'labtest_integrated_numeric.csv',\n",
    "    raw_meta=folder_data / 'ald_metadata_rawfiles.csv')\n",
    "fnames = pimmslearn.nb.Config.from_dict(fnames)  # could be handeled kwargs as in normal dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f38385c7-e42a-40ca-8d37-b7c9c94f4a2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0d19ef8-f566-40dd-ba6c-d47c68d501c5",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43a9b547-a31a-46e8-9a2a-803b936d4738",
   "metadata": {},
   "outputs": [],
   "source": [
    "VAR_PEP = 'PEP.Quantity'\n",
    "VAR_PG = 'PG.Quantity'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0090ae47-f574-4242-a7d3-98201a70d6fa",
   "metadata": {},
   "source": [
    "# Meta data\n",
    "\n",
    "- sample annotation (to select correct samples)\n",
    "- clinical data\n",
    "- meta data from raw files (MS machine recorded meta data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f52877-df54-4c0d-baef-9dcb7c941b4d",
   "metadata": {},
   "source": [
    "## Sample annotations\n",
    "\n",
    "- `Groups`: more detailed (contains sub-batch information)\n",
    "- `Group2`: used to separate samples into cohorts for study\n",
    "- `Sample type`: There are liver biopsy samples measured -> select only Plasma samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c660fa74-9cf1-4aff-b8c1-9e55d5cc38a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "annotations = pd.read_csv(fnames.annotations, index_col='Sample ID')\n",
    "annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9852e507-f636-477b-b539-6e02b6a23b47",
   "metadata": {},
   "outputs": [],
   "source": [
    "annotations['Participant ID'].value_counts().value_counts()  # some only have a blood sample, some both"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1524dead-b5bf-4ee7-acce-6a2dfad890e5",
   "metadata": {},
   "source": [
    "### Select ALD subcohort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad24b48e-9a3f-4c51-9d73-b2a72ed35145",
   "metadata": {},
   "outputs": [],
   "source": [
    "groups = ['ALD']  # 'ALD-validation', 'HP'\n",
    "\n",
    "# annotations.Groups.value_counts()\n",
    "annotations.Group2.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d64210f8-4a44-48fb-8dce-2a2759484e4f",
   "metadata": {},
   "source": [
    "### Select plasma samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11ec53b9-0164-4489-a3fd-38124fc63be8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sel_plasma_samples = (annotations.Group2.isin(['ALD'])) & (annotations['Sample type'] == 'Plasma')\n",
    "sel_plasma_samples = sel_plasma_samples.loc[sel_plasma_samples].index\n",
    "annotations.loc[sel_plasma_samples].describe(include=['object', 'string'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f082cc8-47dd-4b4e-b0bb-60df628fc158",
   "metadata": {},
   "source": [
    "### Select liver samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92633dab-d344-4704-8b4b-658f8c0fe0cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "groups = ['ALD']  # 'ALD-validation', 'HP'\n",
    "sel_liver_samples = (annotations.Group2.isin(['ALD'])) & (annotations['Sample type'] == 'Liver')\n",
    "sel_liver_samples = sel_liver_samples.loc[sel_liver_samples].index\n",
    "annotations.loc[sel_liver_samples].describe(include=['object', 'string'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "952e9a2a-c0fa-426f-a264-f16fd62e1a14",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_qc_plasma = annotations.Group2[annotations.Group2 == 'QC'].index\n",
    "with (folder_data_out / 'qc_samples.yaml').open('w') as f:\n",
    "    yaml.safe_dump(idx_qc_plasma.to_list(), f)\n",
    "idx_qc_plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1b07f3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_qc_liver = annotations.Group2[annotations.Group2 == 'QC_liver'].index\n",
    "with (folder_data_out / 'qc_samples.yaml').open('w') as f:\n",
    "    yaml.safe_dump(idx_qc_liver.to_list(), f)\n",
    "idx_qc_liver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c33b7621-e672-42f5-b3e5-472d68e437ad",
   "metadata": {},
   "source": [
    "## Clinical data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a67a35-f1fd-485f-91ad-2c66563d52dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic = pd.read_csv(fnames.clinic, index_col=0)\n",
    "clinic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6721cb5-b27e-4513-a596-c9be269262d3",
   "metadata": {},
   "source": [
    "- `idx_overlap_plasma`:  Will be used to select samples with data across datasets available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45675ffe-50bd-47fd-8562-7b973fd804ab",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "print('Missing labels: ', sel_plasma_samples.difference(clinic.index))\n",
    "idx_overlap_plasma = clinic.index.intersection(sel_plasma_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e85ecd8e-f1dd-4b3a-81e8-0e28b1deb405",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic.loc[idx_overlap_plasma]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0be1d40a",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic['abstinent_num'] = (clinic[\"currentalc\"] == 0.00).astype(int)\n",
    "clinic[['abstinent_num', 'currentalc']].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c59939eb-7318-4dbd-b00a-40ba4b273f0e",
   "metadata": {},
   "source": [
    "Kleiner score of 0.5 was assigned as value of 0-1 without biopsy. Is set to NA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "796f8882-3a61-41f2-b5a8-da21861a5b04",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic[\"kleiner\"] = clinic[\"kleiner\"].replace({-1: np.nan, 0.5: np.nan})\n",
    "clinic[\"kleiner\"].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1f7a4d5-8ab6-4b07-83cf-27181272ea38",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinic.loc[idx_overlap_plasma].to_csv(folder_data_out / 'ald_metadata_cli.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47d1d427-8d36-44be-8fcb-1889a4aad157",
   "metadata": {},
   "source": [
    "## Rawfile information\n",
    "\n",
    "- [ ] liver samples are currently missing\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fffc83c3-1ee0-49eb-ab91-dba16ef5394c",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_meta = pd.read_csv(fnames.raw_meta, header=[0, 1], index_col=0)\n",
    "raw_meta.index.name = \"Sample ID (long)\"\n",
    "raw_meta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18bfc336-5b95-41b4-8019-b806a1fd0db3",
   "metadata": {},
   "source": [
    "Measurements are super homogenous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aa3aef4-cb6f-4795-9ca1-6357b03f0375",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_meta.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c805b16d-e65b-470e-af34-aa8d7c39e9d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = raw_meta.index.to_series()\n",
    "idx = idx.str.extract(r'(Plate[\\d]_[A-H]\\d*)').squeeze()\n",
    "idx.name = 'Sample ID'\n",
    "idx.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55e66513-cd53-4bdd-940f-1bdfd87e4da2",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_meta = raw_meta.set_index(idx)\n",
    "raw_meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d452fc09-155b-4b48-9498-7bde9658fe86",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_meta_rawfiles_columns = raw_meta.columns  # needs to go to Config which is not overwriteable by attribute selection\n",
    "meta_raw_names = raw_meta.columns.droplevel()\n",
    "assert meta_raw_names.is_unique\n",
    "meta_raw_names.name = None\n",
    "raw_meta.columns = meta_raw_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1996dd1-8dca-4bd4-b7c5-a3b41b20ddfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_meta.loc[['Plate6_F2']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80878aca-78b3-4533-b338-a3fa935a7763",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Missing metadata in set of selected labels: \", idx_overlap_plasma.difference(raw_meta.index))\n",
    "idx_overlap_plasma = idx_overlap_plasma.intersection(raw_meta.index)  # proteomics data has to be part of metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c09d2356-99bf-43f4-ae9b-831f0d011c92",
   "metadata": {},
   "source": [
    "Still save all metadata which is there, but subselect data samples accordingly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bedb3ea-454a-41eb-9b5b-576f39185cd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_meta.to_csv(folder_data_out / 'raw_meta.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baa62399-a23a-44e8-8479-1a8c2aa587df",
   "metadata": {},
   "source": [
    "# Plasma samples\n",
    "\n",
    "- load samples\n",
    "- select based on `sel_plasma_samples`, `raw_meta`  (inclusion in clinical cohort was done before)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36f9abaf-4e9c-48f0-ab7c-ebb12c06dbad",
   "metadata": {},
   "source": [
    "## Missing samples\n",
    "\n",
    "From the above we can note that there is\n",
    "- no clinical data for `Plate6_F2`\n",
    "- no metadata for `Plate2_C1`: re-measured sample which looks fine, but fails with error `\"Unable to access the RAW file using the native Thermo library\"`\n",
    "\n",
    "> see section below"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfe1c458-dc61-4890-b430-6efa7eb89e72",
   "metadata": {},
   "source": [
    "## (Aggregated) Peptide Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5909039c-1dec-4aaf-a889-d43d33e8a324",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_table(fnames.plasma_aggPeptides, low_memory=False)\n",
    "N_FRIST_META = 8\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e15d6b3-2735-4f82-9eff-09097bbcc72e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1faea089",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[:, :N_FRIST_META].describe(include='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd7b6bbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = df.iloc[:, N_FRIST_META:].columns.to_series().apply(lambda s: tuple(s.split('.')[-2:]))\n",
    "column_types.describe()  # .apply(lambda l: l[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96399cdc-beee-45ac-8f65-2bd7ca9d9f9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = ['.'.join(x for x in tup) for tup in list(column_types.unique())]\n",
    "column_types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aa3fccb-57f7-47fc-a6aa-9028972cacee",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(list(df.columns[:N_FRIST_META])).sort_index(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d12c9ec4-66e7-4df5-bfa3-62e570c0ebf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:, df.columns.str.contains(VAR_PEP)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5203a536-3838-44a0-bd18-b9d3b8a645af",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[:20, :6]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e471878b-6198-49d7-848f-cb2f62c56c95",
   "metadata": {},
   "source": [
    "create new multiindex from column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faca2822",
   "metadata": {},
   "outputs": [],
   "source": [
    "sep = '.raw.'\n",
    "df.columns = pd.MultiIndex.from_tuples(df.columns.str.split().str[1].str.split(\n",
    "    sep).to_series().apply(tuple), names=['Sample ID', 'vars'])\n",
    "df = df.stack(0)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4698cfd0-59f0-44cd-aefa-1a5febc84a1e",
   "metadata": {},
   "source": [
    "### Index meta data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b06426c-9f19-4602-9176-55778df36ead",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "meta = df.index.to_frame().reset_index(drop=True)\n",
    "meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0964a2d2-2969-4754-9b74-f5a5b7a4d60c",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta.describe(include='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cb30a12",
   "metadata": {},
   "outputs": [],
   "source": [
    "id_mappings = [\"PEP.StrippedSequence\", \"PG.ProteinAccessions\", \"PG.Genes\"]\n",
    "id_mappings = meta[id_mappings].drop_duplicates()\n",
    "id_mappings.to_csv(folder_data_out / 'ald_plasma_aggPeptides_id_mappings.csv')\n",
    "id_mappings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d28918f8-32b7-46af-8545-1f998bb3efc3",
   "metadata": {},
   "source": [
    "### Select aggregated peptide level data\n",
    "\n",
    "taken from [Spectronaut manuel](https://biognosys.com/resources/spectronaut-manual/)\n",
    "\n",
    "feature | description\n",
    "--- | ---\n",
    "PEP.IsProteinGroupSpecific | True or False. Tells you whether the peptide only belongs to one Protein Group.\n",
    "PEP.StrippedSequence | -\n",
    "PEP.IsProteotypic |  -\n",
    "PEP.PeptidePosition | -\n",
    "PG.Cscore | -\n",
    "PG.ProteinAccessions | -\n",
    "PG.Genes | -\n",
    "PEP.Quantity | The quantitative value for that peptide as defined in the settings.\n",
    "EG.PrecursorId | Unique Id for the precursor: [modified sequence] plus [charge]\n",
    "EG.Qvalue | The q-value (FDR) of the EG.\n",
    "EG.TotalQuantity (Settings) | The quantitative value for that EG as defined in the settings.\n",
    "\n",
    "> Headers related to Peptides (PEP) as defined in the settings. Many headers related to Peptides are self-explanatory. Here are the most relevant and some which are not too obvious.\n",
    "\n",
    "> Headers related to Peptides (PEP) as defined in the settings. Many headers related to Peptides are self-explanatory. Here are the most relevant and some which are not too obvious.\n",
    "\n",
    "After discussing with Lili, `PEP.Quantity` is the fitting entity for\n",
    "each unique aggregated Peptide. Duplicated entries are just to drop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbdb7430-0acd-43dc-a614-4ac3541a6422",
   "metadata": {},
   "outputs": [],
   "source": [
    "sel_cols = ['Sample ID', 'PEP.StrippedSequence', 'PEP.Quantity']  # selected quantity in last position\n",
    "df = df.reset_index()[sel_cols].drop_duplicates().set_index(sel_cols[:2])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ee74a73",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.squeeze().dropna().astype(float).unstack()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ebe72f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.dropna(how='all', axis=1)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f1bdeeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = df.index.to_series()\n",
    "idx = idx.str.extract(r'(Plate[\\d]_[A-H]\\d*)').squeeze()\n",
    "idx.name = 'Sample ID'\n",
    "idx.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b823acaf-2610-4b0a-91d8-2d6dd6ff4182",
   "metadata": {},
   "source": [
    "- rawfile metadata -> keep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eb9f6d4-8f0b-40d5-9913-67e9119124a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(idx)\n",
    "df_qc = df.loc[idx_qc_plasma].copy()\n",
    "df = df.loc[idx_overlap_plasma]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5ddb733",
   "metadata": {},
   "outputs": [],
   "source": [
    "# des_data = df.describe() # too slow\n",
    "des_data = df.isna().sum().to_frame('count').T\n",
    "des_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8591f292-e983-4fce-90bb-f6dc75c5b482",
   "metadata": {},
   "source": [
    "### Check for metadata from rawfile overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dddafbb-edd7-4ef0-9787-3120b24d7f79",
   "metadata": {},
   "source": [
    "For one raw file no metadata could be extracted (`ERROR: Unable to\n",
    "access the RAW file using the native Thermo library.`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8a09c6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_diff = df.index.difference(raw_meta.index)\n",
    "annotations.loc[idx_diff]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9037f332",
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs = {'xlabel': 'peptide number ordered by completeness',\n",
    "          'ylabel': 'peptide was found in # samples',\n",
    "          'title': 'peptide measurement distribution'}\n",
    "\n",
    "ax = pimmslearn.plotting.plot_counts(des_data.T.sort_values(by='count', ascending=False).reset_index(\n",
    "), feat_col_name='count', feature_name='Aggregated peptides', n_samples=len(df), ax=None, **kwargs)\n",
    "\n",
    "fig = ax.get_figure()\n",
    "fig.tight_layout()\n",
    "pimmslearn.savefig(fig, name='data_aggPeptides_completness', folder=folder_run)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9bf0653-6437-47c3-b41a-47c3eb41b6b0",
   "metadata": {},
   "source": [
    "### Select features which are present in at least 25% of the samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee502baa",
   "metadata": {},
   "outputs": [],
   "source": [
    "PROP_FEAT_OVER_SAMPLES = .25\n",
    "prop = des_data.loc['count'] / len(df)\n",
    "selected = prop >= PROP_FEAT_OVER_SAMPLES\n",
    "selected.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c25e56b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.loc[:, selected]\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e78b6f3e-a097-4891-b7f2-2f9c185f226f",
   "metadata": {},
   "source": [
    "Dump selected data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3b25708-a88b-4c1d-9271-016614160dc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames.sel_plasma_aggPeptides = folder_data_out / 'ald_plasma_aggPeptides.pkl'\n",
    "df.to_pickle(fnames.sel_plasma_aggPeptides)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83c4e473-6d94-46e1-954c-ddadfbc2b188",
   "metadata": {},
   "source": [
    "Dump QC sample data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d468a3e-f186-464d-9fcb-21b7aa455be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_qc = df_qc.loc[:, selected]\n",
    "fnames.qc_plasma_aggPeptides = folder_data_out / 'qc_plasma_aggPeptides.pkl'\n",
    "df_qc.to_pickle(fnames.qc_plasma_aggPeptides)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ebd784a-9c35-41f3-beb0-711f8475a742",
   "metadata": {},
   "source": [
    "## Protein Group data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5ef1acb-1cb4-4e08-bad9-a44aef9c2a58",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(fnames.plasma_proteinGroups, low_memory=False)\n",
    "idx_cols = ['PG.ProteinAccessions', 'PG.Genes']\n",
    "N_FRIST_META = 3\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cf7a8e2-0540-49b4-97c3-d9c8efb59f8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta = df.iloc[:, :N_FRIST_META]\n",
    "meta.describe(include='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48e14dd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "id_mappings = [\"PG.ProteinAccessions\", \"PG.Genes\"]\n",
    "id_mappings = meta[id_mappings].drop_duplicates()\n",
    "id_mappings.to_csv(folder_data_out / 'ald_plasma_proteinGroups_id_mappings.csv', index=False)\n",
    "id_mappings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b811c28-d282-4c53-93a0-febef823e537",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = df.iloc[:, N_FRIST_META:].columns.to_series().apply(lambda s: tuple(s.split('.')[-2:]))\n",
    "column_types.describe()  # .apply(lambda l: l[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e9611d6-575e-42ac-ac2f-650532b86b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = ['.'.join(x for x in tup) for tup in list(column_types.unique())]\n",
    "column_types  # 'PG.Quantity' expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00f2c254",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(list(df.columns[:N_FRIST_META])).sort_index(axis=1)\n",
    "df.loc[:, df.columns.str.contains(VAR_PG)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97272d7d-3c81-441c-b45e-d79cd858c96d",
   "metadata": {},
   "source": [
    "Drop index columns which are not selected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ad74f5e-2555-431b-bb88-ac65891f816f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_idx_to_drop(df: pd.DataFrame, idx_to_keep: list):\n",
    "    to_drop = [x for x in df.index.names if x not in idx_to_keep]\n",
    "    logger.info(\"Columnns to drop: {}\".format(\",\".join((str(x) for x in to_drop))))\n",
    "    return to_drop\n",
    "\n",
    "\n",
    "to_drop = find_idx_to_drop(df, idx_cols)\n",
    "df = df.reset_index(level=to_drop, drop=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "257f91cc-a833-4601-aa0d-1c4b180c0067",
   "metadata": {},
   "source": [
    "extract long sample name (highly specific to task)\n",
    "- whitespace split, taking last position of column name\n",
    "- `sep` splits `Sample ID` from `vars`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75cbb2f8-7fde-421d-9b3d-d7ce1dfc5052",
   "metadata": {},
   "outputs": [],
   "source": [
    "sep = '.raw.'\n",
    "# sep = '.htrms.'\n",
    "df.columns = pd.MultiIndex.from_tuples(df.columns.str.split().str[-1].str.split(\n",
    "    sep).to_series().apply(tuple), names=['Sample ID', 'vars'])\n",
    "df = df.stack(0)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cbaf13b-437c-4e93-8712-4cc3afb0497e",
   "metadata": {},
   "source": [
    "### Select Protein Group data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f5c011",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[[VAR_PG]]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "041cb10d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = df['PG.Quantity'] == 'Filtered'\n",
    "print(\"No. of Filtered entries: \", mask.sum())\n",
    "df = df.loc[~mask]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8042ccc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cdfec5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.squeeze().dropna().astype(float).unstack()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b291f117",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_non_unique = df.index.to_frame()[\"PG.Genes\"].value_counts() > 1\n",
    "gene_non_unique = gene_non_unique[gene_non_unique].index\n",
    "gene_non_unique"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea188e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[pd.IndexSlice[:, gene_non_unique], :].T.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7e1196f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.T\n",
    "\n",
    "idx = df.index.to_series()\n",
    "idx = idx.str.extract(r'(Plate[\\d]_[A-H]\\d*)').squeeze()\n",
    "idx.name = 'Sample ID'\n",
    "idx.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9878be79",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(idx)\n",
    "df_qc = df.loc[idx_qc_plasma].copy()\n",
    "df = df.loc[idx_overlap_plasma]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7aec3aaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.dropna(how='all', axis=0)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57a6f3c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "des_data = df.describe()\n",
    "des_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ffb7e87-befc-4a34-b96e-4609dc23b8fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_feat = des_data.loc[\"count\"].droplevel(-1).rename('freq')\n",
    "freq_feat.to_csv(folder_data_out / 'freq_ald_plasma_proteinGroups.csv')\n",
    "# pd.read_csv(folder_data_out / 'freq_ald_plasma_proteinGroups.csv', index_col=0)\n",
    "freq_feat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f31ff36-c456-427d-b2de-03f43bb51f16",
   "metadata": {},
   "source": [
    "### Check for metadata from rawfile overlap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abb0778c-e611-476c-9272-78fd4ee35afc",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_diff = df.index.difference(raw_meta.index)\n",
    "annotations.loc[idx_diff]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8967480e-9657-4c9d-b42a-fb88e03c61b2",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "kwargs = {'xlabel': 'protein group number ordered by completeness',\n",
    "          'ylabel': 'peptide was found in # samples',\n",
    "          'title': 'protein group measurement distribution'}\n",
    "\n",
    "ax = pimmslearn.plotting.plot_counts(des_data.T.sort_values(by='count', ascending=False).reset_index(\n",
    "), feat_col_name='count', n_samples=len(df), ax=None, min_feat_prop=.0, **kwargs)\n",
    "\n",
    "fig = ax.get_figure()\n",
    "fig.tight_layout()\n",
    "pimmslearn.savefig(fig, name='data_proteinGroups_completness', folder=folder_run)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1155349a",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "Save unfiltered data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef131086",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "fname = folder_data_out / 'ald_plasma_proteinGroups_unfiltered.pkl'\n",
    "df.to_pickle(fname)\n",
    "fname"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b64a1eb5-32ad-4b75-8dc9-0de6c873ff3d",
   "metadata": {},
   "source": [
    "### Select features which are present in at least 25% of the samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05f024ee-54bc-452e-beff-9518a82471c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "PROP_FEAT_OVER_SAMPLES = .25\n",
    "prop = des_data.loc['count'] / len(df)\n",
    "selected = prop >= PROP_FEAT_OVER_SAMPLES\n",
    "selected.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36bdad90-eae3-48d2-bbc8-87c65943cab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.loc[:, selected]\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df435a1-8e7e-4888-8bae-473748a1d9f2",
   "metadata": {},
   "source": [
    "Check for non unique genes after dropping uncommon protein groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30121f05",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_non_unique = df.columns.to_frame()[\"PG.Genes\"].value_counts() > 1\n",
    "gene_non_unique = gene_non_unique[gene_non_unique].index\n",
    "gene_non_unique"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bd35030-567c-414f-928a-6f643aa9fb33",
   "metadata": {},
   "source": [
    "- less often found -> less intensity on average and on maximum\n",
    "\n",
    "- [ ] decided if protein group should be subselected\n",
    "- alternative selection: per sample, select protein group with highest intensity per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "132cf838",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.T.loc[pd.IndexSlice[:, gene_non_unique], :].T.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2f3adaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.droplevel(1, axis=1)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "073ffa7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_pickle(folder_data_out / 'ald_plasma_proteinGroups.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df3ad9ff-4f2a-4aca-958b-bd963f0d58a1",
   "metadata": {},
   "source": [
    "Dump QC sample data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "737aeda4-be8b-4513-a636-949fedbe1c61",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_qc = df_qc.loc[:, selected].droplevel(1, axis=1)\n",
    "fnames.qc_plasma_proteinGroups = folder_data_out / 'qc_plasma_proteinGroups.pkl'\n",
    "df_qc.to_pickle(fnames.qc_plasma_proteinGroups)\n",
    "df_qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "857bdb56-e185-4cd9-82a0-f25dca0dcc59",
   "metadata": {},
   "source": [
    "# Liver samples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e316b53-64d6-40b5-b6f3-9c01d6477748",
   "metadata": {},
   "source": [
    "## Peptides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77c1829d",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_cols = ['PG.ProteinAccessions', 'PG.Genes', 'Sample ID']\n",
    "N_FRIST_META = 8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0641f1a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(fnames.liver_aggPeptides, low_memory=False)\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a39b468c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a91135f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[:, :N_FRIST_META].describe(include='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27b62e69",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = df.iloc[:, N_FRIST_META:].columns.to_series().apply(lambda s: tuple(s.split('.')[-2:]))\n",
    "column_types.describe()  # .apply(lambda l: l[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e5e1416",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = ['.'.join(x for x in tup) for tup in list(column_types.unique())]\n",
    "column_types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a2846a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(list(df.columns[:N_FRIST_META])).sort_index(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7227cd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:, df.columns.str.contains(VAR_PEP)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "603f3cc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b710c37d-2d57-44ec-a1e3-00b02106d412",
   "metadata": {},
   "source": [
    "create new multiindex from column (see examples above)\n",
    "\n",
    "- split on whitespace\n",
    "- select string at first position\n",
    "- split on `sep`, keep both"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d14f3917",
   "metadata": {},
   "outputs": [],
   "source": [
    "sep = '.htrms.'\n",
    "df.columns = pd.MultiIndex.from_tuples(df.columns.str.split().str[1].str.split(\n",
    "    sep).to_series().apply(tuple), names=['Sample ID', 'vars'])\n",
    "df = df.stack(0)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e46a46d9-489f-4d01-9338-00aaacfc0e1a",
   "metadata": {},
   "source": [
    "### Index meta data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d157da80",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "meta = df.index.to_frame().reset_index(drop=True)\n",
    "meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c1e81d8",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "id_mappings = [\"PEP.StrippedSequence\", \"PG.ProteinAccessions\", \"PG.Genes\"]\n",
    "id_mappings = meta[id_mappings].drop_duplicates()\n",
    "id_mappings.to_csv(folder_data_out / 'ald_liver_aggPeptides_id_mappings.csv')\n",
    "id_mappings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b981495b-9305-4df2-9aea-3cd5bfe01eee",
   "metadata": {},
   "source": [
    "### Select aggregated peptide level data\n",
    "\n",
    "taken from [Spectronaut manuel](https://biognosys.com/resources/spectronaut-manual/)\n",
    "\n",
    "feature | description\n",
    "--- | ---\n",
    "PEP.IsProteinGroupSpecific | True or False. Tells you whether the peptide only belongs to one Protein Group.\n",
    "PEP.StrippedSequence | -\n",
    "PEP.IsProteotypic |  -\n",
    "PEP.PeptidePosition | -\n",
    "PG.Cscore | -\n",
    "PG.ProteinAccessions | -\n",
    "PG.Genes | -\n",
    "PEP.Quantity | The quantitative value for that peptide as defined in the settings.\n",
    "EG.PrecursorId | Unique Id for the precursor: [modified sequence] plus [charge]\n",
    "EG.Qvalue | The q-value (FDR) of the EG.\n",
    "EG.TotalQuantity (Settings) | The quantitative value for that EG as defined in the settings.\n",
    "\n",
    "> Headers related to Peptides (PEP) as defined in the settings. Many headers related to Peptides are self-explanatory. Here are the most relevant and some which are not too obvious.\n",
    "\n",
    "> Headers related to Peptides (PEP) as defined in the settings. Many headers related to Peptides are self-explanatory. Here are the most relevant and some which are not too obvious.\n",
    "\n",
    "After discussing with Lili, `PEP.Quantity` is the fitting entity for\n",
    "each unique aggregated Peptide. Duplicated entries are just to drop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41eaa101",
   "metadata": {},
   "outputs": [],
   "source": [
    "sel_cols = ['Sample ID', 'PEP.StrippedSequence', VAR_PEP]  # selected quantity in last position\n",
    "df = df.reset_index()[sel_cols].drop_duplicates().set_index(sel_cols[:2]).squeeze()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "365fbb69",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = df != 'Filtered'\n",
    "df = df.loc[mask].astype(float)\n",
    "df.sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fe42ef8-a6ec-4a72-bc24-5c218370c5f4",
   "metadata": {},
   "source": [
    "Select entry with maximum intensity of `duplicated entries`\n",
    "\n",
    "> change of variable and many duplicates -> could be PSM table? (close to evidence?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9047d416-01b9-4faf-ad77-d48d782cab1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_idx_duplicated = df.index.duplicated(False)\n",
    "df.loc[mask_idx_duplicated].sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25fa85cb-f408-4f49-b360-f7c75b5a06c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pimmslearn.pandas.select_max_by(df=df.reset_index(),\n",
    "                               grouping_columns=sel_cols[:-1],\n",
    "                               selection_column=sel_cols[-1]).set_index(sel_cols[:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb0a4937-525d-4338-ad2f-f9dbc74e1010",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert df.index.duplicated(False).sum() == 0, \"Still missing values\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "701a04c0-ab60-4ab4-afc6-12920807ead0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.unstack()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb39ab1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = df.index.to_series()\n",
    "idx = idx.str.extract(r'(PlateS[\\d]_[A-H]\\d*)').squeeze()\n",
    "idx.name = 'Sample ID'\n",
    "idx.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529fa0e7-7ad4-4c72-91b1-d37587835ce5",
   "metadata": {},
   "source": [
    "- rawfile metadata -> keep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f987b598",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(idx)\n",
    "df_qc = df.loc[idx_qc_liver]\n",
    "df = df.loc[sel_liver_samples]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f92a16e1-51d8-4b55-977a-3627c9abd35c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.dropna(how='all', axis=1)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0671ad9-6b37-4cbc-a636-bab469253c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# des_data = df.describe() unnecessary computation which take too long\n",
    "des_data = df.isna().sum().to_frame('count').T\n",
    "des_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a12c83be-dc97-4604-bc88-2ea97b058e19",
   "metadata": {},
   "source": [
    "### Check for metadata from rawfile overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44616770-fcc2-4a97-86f4-e0eadc98bb7a",
   "metadata": {},
   "source": [
    "For one raw file no metadata could be extracted (`ERROR: Unable to\n",
    "access the RAW file using the native Thermo library.`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfc5e770",
   "metadata": {},
   "outputs": [],
   "source": [
    "# idx_diff = df.index.difference(raw_meta.index)\n",
    "# annotations.loc[idx_diff]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2ca9fcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs = {'xlabel': 'peptide number ordered by completeness',\n",
    "          'ylabel': 'peptide was found in # samples',\n",
    "          'title': 'peptide measurement distribution'}\n",
    "\n",
    "ax = pimmslearn.plotting.plot_counts(des_data.T.sort_values(by='count', ascending=False).reset_index(\n",
    "), feat_col_name='count', feature_name='Aggregated peptides', n_samples=len(df), ax=None, **kwargs)\n",
    "\n",
    "fig = ax.get_figure()\n",
    "fig.tight_layout()\n",
    "pimmslearn.savefig(fig, name='data_liver_aggPeptides_completness', folder=folder_run)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "920ee2e1-d99c-4203-8c84-e5a02b73770f",
   "metadata": {},
   "source": [
    "### Select features which are present in at least 25% of the samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddf9c82e",
   "metadata": {},
   "outputs": [],
   "source": [
    "PROP_FEAT_OVER_SAMPLES = .25\n",
    "prop = des_data.loc['count'] / len(df)\n",
    "selected = prop >= PROP_FEAT_OVER_SAMPLES\n",
    "selected.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "000e2afb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.loc[:, selected]\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1467b886-bf4d-4f31-976e-cd8edc3d3b41",
   "metadata": {},
   "source": [
    "Dump selected data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d79ce44c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames.sel_liver_aggPeptides = folder_data_out / 'ald_liver_aggPeptides.pkl'\n",
    "df.to_pickle(fnames.sel_liver_aggPeptides)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ba844aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames.qc_liver_aggPeptides = folder_data_out / 'qc_liver_aggPeptides.pkl'\n",
    "df_qc.to_pickle(fnames.qc_liver_aggPeptides)\n",
    "df_qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cd3f6a6-36b8-49b2-9f21-9298446e68ff",
   "metadata": {},
   "source": [
    "## Protein Groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92494f1c-231b-4fb6-bf4d-5096f2ecc07b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(fnames.liver_proteinGroups, low_memory=False)\n",
    "idx_cols = ['PG.ProteinAccessions', 'PG.Genes']\n",
    "N_FRIST_META = 5\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5b83ac5-44da-4df8-b08e-5cfd40b88822",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# find_idx_to_drop(df, idx_cols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01ccd4e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta = df.iloc[:, :N_FRIST_META]\n",
    "meta.describe(include='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01ed7326",
   "metadata": {},
   "outputs": [],
   "source": [
    "id_mappings = [\"PG.ProteinAccessions\", \"PG.Genes\"]\n",
    "id_mappings = meta[id_mappings].drop_duplicates()\n",
    "id_mappings.to_csv(folder_data_out / 'ald_liver_proteinGroups_id_mappings.csv')\n",
    "id_mappings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43470ad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = df.iloc[:, N_FRIST_META:].columns.to_series().apply(lambda s: tuple(s.split('.')[-2:]))\n",
    "column_types.describe()  # .apply(lambda l: l[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ff5e4c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_types = ['.'.join(x for x in tup) for tup in list(column_types.unique())]\n",
    "column_types  # 'PG.Quantity' expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4154c47",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "df = df.set_index(list(df.columns[:N_FRIST_META])).sort_index(axis=1)\n",
    "df.loc[:, df.columns.str.contains(VAR_PG)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4e32e4",
   "metadata": {},
   "source": [
    "Drop index columns which are not selected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f56d6c6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to_drop = find_idx_to_drop(df, idx_cols)\n",
    "# df = df.reset_index(level=to_drop, drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a372506",
   "metadata": {},
   "source": [
    "extract long sample name (highly specific to task)\n",
    "- whitespace split, taking last position of column name\n",
    "- `sep` splits `Sample ID` from `vars`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ba2128b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sep = '.htrms.'\n",
    "df.columns = pd.MultiIndex.from_tuples(df.columns.str.split().str[-1].str.split(\n",
    "    sep).to_series().apply(tuple), names=['Sample ID', 'vars'])\n",
    "df = df.stack(0)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bac64d1b",
   "metadata": {},
   "source": [
    "### Select Protein Group data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c70ff1fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[[VAR_PG]]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "711002a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = df['PG.Quantity'] == 'Filtered'\n",
    "print(\"No. of Filtered entries: \", mask.sum())\n",
    "df = df.loc[~mask]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "440ea167-d84e-4c8a-89a7-2ad47665c6f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sel_cols = ['PG.ProteinAccessions', 'PG.Genes', 'Sample ID', VAR_PG]  # last one gives quantity\n",
    "df = df.reset_index()[sel_cols].drop_duplicates().set_index(sel_cols[:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "181137fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d38ffb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.squeeze().dropna().astype(float).unstack()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fb9660a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_non_unique = df.index.to_frame()[\"PG.Genes\"].value_counts() > 1\n",
    "gene_non_unique = gene_non_unique[gene_non_unique].index\n",
    "gene_non_unique"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b70499c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[pd.IndexSlice[:, gene_non_unique], :].T.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94c7d042",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.T\n",
    "\n",
    "idx = df.index.to_series()\n",
    "idx = idx.str.extract(r'(PlateS[\\d]_[A-H]\\d*)').squeeze()\n",
    "idx.name = 'Sample ID'\n",
    "idx.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71c45282",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(idx)\n",
    "df_qc = df.loc[idx_qc_liver]\n",
    "df = df.loc[sel_liver_samples]\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a1c2dbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.dropna(how='all', axis=1)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85ba0421",
   "metadata": {},
   "outputs": [],
   "source": [
    "des_data = df.describe()\n",
    "des_data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f15400ca",
   "metadata": {},
   "source": [
    "### Check for metadata from rawfile overlap\n",
    "- no raw data yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae6f31d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_diff = df.index.difference(raw_meta.index)\n",
    "annotations.loc[idx_diff]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c66f5976",
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs = {'xlabel': 'protein group number ordered by completeness',\n",
    "          'ylabel': 'peptide was found in # samples',\n",
    "          'title': 'protein group measurement distribution'}\n",
    "\n",
    "ax = pimmslearn.plotting.plot_counts(des_data.T.sort_values(by='count', ascending=False).reset_index(\n",
    "), feat_col_name='count', n_samples=len(df), ax=None, **kwargs)\n",
    "\n",
    "fig = ax.get_figure()\n",
    "fig.tight_layout()\n",
    "fnames.fig_liver_pg_completness = folder_run / 'data_liver_proteinGroups_completness'\n",
    "pimmslearn.savefig(fig, name=fnames.fig_liver_pg_completness)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "205fa9fd",
   "metadata": {},
   "source": [
    "### Select features which are present in at least 25% of the samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9e03b6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "PROP_FEAT_OVER_SAMPLES = .25\n",
    "prop = des_data.loc['count'] / len(df)\n",
    "selected = prop >= PROP_FEAT_OVER_SAMPLES\n",
    "selected.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d76569d",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.loc[:, selected]\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c46c0973",
   "metadata": {},
   "source": [
    "Check for non unique genes after dropping uncommon protein groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9191dd2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_non_unique = df.columns.to_frame()[\"PG.Genes\"].value_counts() > 1\n",
    "gene_non_unique = gene_non_unique[gene_non_unique].index\n",
    "gene_non_unique"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf8be578",
   "metadata": {},
   "source": [
    "- less often found -> less intensity on average and on maximum\n",
    "\n",
    "- [ ] decided if protein group should be subselected\n",
    "- alternative selection: per sample, select protein group with highest intensity per sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b1eff86",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.T.loc[pd.IndexSlice[:, gene_non_unique], :].T.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8971212",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.droplevel(1, axis=1)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7a1067a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames.sel_liver_proteinGroups = folder_data_out / 'ald_liver_proteinGroups.pkl'\n",
    "df.to_pickle(fnames.sel_liver_proteinGroups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51ee65d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames.qc_liver_proteinGroups = folder_data_out / 'qc_liver_proteinGroups.pkl'\n",
    "df_qc.to_pickle(fnames.qc_liver_proteinGroups)\n",
    "df_qc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3a38073-d669-4859-9308-e0f7bc16f0a2",
   "metadata": {},
   "source": [
    "# All file names\n",
    "- inputs\n",
    "- output data (`pkl`)\n",
    "- figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3e2b4eb-bf65-46be-8568-c8f67accda0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9013cc81-684b-4425-8775-e409498d47e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfg_dump = 'config/ald_data.yaml'\n",
    "fnames.dump(cfg_dump)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}