ALD-PA/ALD_PA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Project and this notebook\n",
    "- A paired liver biopsy and plasma proteomics study reveals circulating biomarkers for alcohol-related liver disease\n",
    "- bioRxiv preprint: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- This notebook contains proteomics data anlaysis in the ALD study that are performed in Python scripts. For the machine learning pipeline, please refer to ALD-ML.ipynb in the same repository. Due to restrictions with sensitive data, clinical and proteomics datasets cannot be made public yet. A preprint of the manuscript is available at bioRxiv: In particular, this notebook contains the following analysis:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overview of datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plasma proteome data from 596 study subjects + 7 QC samples\n",
    "- Liver biopsy proteome data from 79 study subjects + 4 QC samples\n",
    "- Clinical data from 595 study subjects \n",
    "- All datasets cannot be made public until the approval for data closure is granted from the Danish Data Protection Agency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overview of analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Experimental design\n",
    "- ID mapping\n",
    "- Clinical data pre-processing\n",
    "- Plasma proteomics data processing\n",
    "- Plasma proteomics data quality \n",
    "- Liver proteomics data processing\n",
    "- Liver proteomics data quality \n",
    "- ANCOVA in the liver proteome\n",
    "- ANCOVA in the plasma proteome\n",
    "- Liver-plasma proteome correlation\n",
    "- Cohort characteristics\n",
    "- Distribution of disease severity in the plasma samples\n",
    "- Distribution of disease severity in the liver samples\n",
    "- Partial correlation in the liver- and plasma proteomes\n",
    "- Global correlation between the liver- and plasma with clinical data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pingouin as pg\n",
    "import scipy.stats as stats\n",
    "from math import log\n",
    "import statsmodels.stats.multitest as multi\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "import seaborn as sns\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def rename_experiment (data, base_label, new_label):\n",
    "    labelfile = pd.read_csv('data/raw/Experiment annotation file_upgrade.csv')\n",
    "    label_dict=dict(zip(labelfile[base_label],labelfile[new_label]))\n",
    "    data=data.rename(columns=label_dict)\n",
    "    return data\n",
    "\n",
    "def extract_datamatrix_DIA(Report):\n",
    "    data = Report.copy()\n",
    "    for i in ['PG.ProteinAccessions', 'PG.Genes']:\n",
    "        data[i]=data[i].str.split(';').str[0]\n",
    "    data=data.rename({'PG.Genes': 'Gene names', 'PG.ProteinAccessions': 'Protein ID'}, axis=1)\n",
    "    data=rename_experiment(data, base_label='File name', new_label='Sample ID')\n",
    "    data=data.replace({'Filtered': np.float64('NaN')})\n",
    "    return data\n",
    "\n",
    "def log_transform_data(df, logbase):\n",
    "    data = df.copy()\n",
    "    if logbase == 2:\n",
    "        data = np.log2(data)\n",
    "    elif logbase == 10:\n",
    "        data = np.log10(data)\n",
    "    else:\n",
    "        data = \"Log base not included\"\n",
    "        \n",
    "    return data\n",
    "\n",
    "def imputation_normal_distribution(df):\n",
    "    data_imputed = df.copy()\n",
    "    for i in data_imputed.loc[:, data_imputed.isnull().any()]:\n",
    "        missing = data_imputed[i].isnull()\n",
    "        std = data_imputed[i].std()\n",
    "        mean = data_imputed[i].mean()\n",
    "        sigma = std*0.3\n",
    "        mu = mean - (std*1.8)\n",
    "        np.random.seed(123)\n",
    "        data_imputed.loc[missing, i] = np.random.normal(mu, sigma, size=len(data_imputed[missing]))\n",
    "        \n",
    "    return data_imputed \n",
    "\n",
    "def tukey_HSD (data, dv, between, fdr = 0.05):\n",
    "    \"\"\"\n",
    "    This is a wrapper of the pingouin.pairwise_tukey for multiple hypothesis test, with fdr control by Benjamini-Hochberg.\n",
    "    Long data format, with protein ID as index.\n",
    "    Pair-wise Tukey-HSD post hoc test, groups should have equal variances and optimally balanced.\n",
    "    When the sample sizes are unequal, this function actually performs the Tukey-Kramer test (which allows for unequal sample sizes).\n",
    "    If the groups have unequal variances, the Games-Howell test is more adequate.\n",
    "    \"\"\"\n",
    "    columns = ['A', 'B', 'mean(A)', 'mean(B)', 'diff', 'SE', 'tail', 'T', 'p-tukey', 'efsize', 'eftype']\n",
    "    scores = pd.DataFrame(columns = columns)\n",
    "    for i in list(set(data.index)):\n",
    "        df_tukeyHSD = data.loc[i]\n",
    "        result = pg.pairwise_tukey(dv = dv, between = between, data = df_tukeyHSD)\n",
    "        result['protein'] = i\n",
    "        scores = scores.append(result)\n",
    "    scores = scores.assign(new_column = lambda x: -np.log10(scores['p-tukey']))\n",
    "    scores = scores.rename({'new_column' : '-Log pvalue'}, axis = 1)\n",
    "    \n",
    "    #FDR correction\n",
    "    np.random.seed(123)\n",
    "    reject, qvalue = multi.fdrcorrection(scores['p-tukey'], alpha=0.05, method='indep')\n",
    "    scores['qvalue'] = qvalue\n",
    "    scores['rejected'] = reject\n",
    "    scores = scores.set_index('protein')\n",
    "    scores['Gene name'] = scores.index.map(IDmapping_UniprotID_to_Genename)\n",
    "    return scores[scores.rejected]\n",
    "\n",
    "def pairwise_correlation (data, columns, fdr = 0.05):\n",
    "    \"\"\"\n",
    "    This is a wrapper of the pingouin.pairwise_corr for multiple hypothesis test, with fdr control by Benjamini-Hochberg.\n",
    "    Long data format, with protein ID as index, pair-wise values as two columns.\n",
    "    \"\"\"\n",
    "    df = data.copy()\n",
    "    columns_new = ['X', 'Y', 'method', 'tail', 'n', 'r', 'CI95%', 'r2', 'adj_r2', 'z','p-unc', 'BF10', 'power']\n",
    "    scores = pd.DataFrame(columns = columns_new)\n",
    "    for i in list(set(data.index)):\n",
    "        df_pairwise_corr = df.loc[i].dropna()\n",
    "        if df_pairwise_corr.shape[0] > 4:\n",
    "            result = pg.pairwise_corr(data = df_pairwise_corr, columns = columns)\n",
    "            result['protein'] = i\n",
    "            scores = scores.append(result)\n",
    "    scores = scores.assign(new_column = lambda x: -np.log10(scores['p-unc']))\n",
    "    scores = scores.rename({'new_column' : '-Log pvalue'}, axis = 1)\n",
    "    \n",
    "    #FDR correction\n",
    "    np.random.seed(123)\n",
    "    reject, qvalue = multi.fdrcorrection(scores['p-unc'], alpha=0.05, method='indep')\n",
    "    scores['qvalue'] = qvalue\n",
    "    scores['rejected'] = reject\n",
    "    scores = scores.set_index('protein')\n",
    "    scores['Gene name'] = scores.index.map(IDmapping_UniprotID_to_Genename)\n",
    "    return scores\n",
    "\n",
    "def convert_to_numeric(data):\n",
    "    df = data.copy()\n",
    "    columns = df.columns\n",
    "    df_new = pd.DataFrame(columns = columns)\n",
    "    for i in columns:\n",
    "        new_values = pd.to_numeric(df[i], errors = 'ignore')\n",
    "        df_new[i] = new_values\n",
    "    return df_new\n",
    "\n",
    "def ancova_pg (data, dv, between, covar, fdr = 0.05):\n",
    "    \"\"\"\n",
    "    Analysis of covariance (ANCOVA) (pg.ancova), long data format, multiple hypothesis testing corrected by Benjamini-Hochberg.\n",
    "    Note that column name of dv shouldn't contain '\\t'\n",
    "    \"data\": should be long data format, with protein ID as index.\n",
    "    \"dv\": Name of column containing the dependant variable.\n",
    "    \"between\": Name of column containing the between factor.\n",
    "    \"covar\": Name(s) of column(s) containing the covariate. \n",
    "    More refer to: https://pingouin-stats.org/generated/pingouin.ancova.html\n",
    "    \"\"\"\n",
    "    columns = ['protein', 'Source', 'SS', 'DF', 'F', 'p-unc']\n",
    "    scores = pd.DataFrame(columns = columns)\n",
    "    for i in list(set(data.index)):\n",
    "        df_ancova = data.loc[i]\n",
    "        ancova = pg.ancova(data = df_ancova, dv = dv, between = between, covar = covar)\n",
    "        num_covar = len(covar)\n",
    "        ancova['protein'] = i\n",
    "        scores = scores.append(ancova)\n",
    "    scores = scores.assign(new_column = lambda x: -np.log10(scores['p-unc']), sort = False)\n",
    "    scores = scores.rename({'new_column' : '-Log pvalue'}, axis = 1)\n",
    "    scores = scores[scores.Source != 'Residual']\n",
    "    \n",
    "    #FDR correction\n",
    "    np.random.seed(123)\n",
    "    reject, qvalue = multi.fdrcorrection(scores['p-unc'], alpha=0.05, method='indep')\n",
    "    scores['qvalue'] = qvalue\n",
    "    scores['rejected'] = reject\n",
    "    return scores\n",
    "\n",
    "def intersection(list1, list2):\n",
    "    new_list = [i for i in list1 if i in list2]\n",
    "    return new_list\n",
    "\n",
    "def partial_corr_pg(data, x_cols, y_col, covar, method='spearman'):\n",
    "    \"\"\"\n",
    "    This is a wrapper for running partial correlation between multiple \"x\"s and single \"y\", correcting for covariates.\n",
    "    The function is based on pingouin.partial_corr. (https://pingouin-stats.org/generated/pingouin.partial_corr.html)\n",
    "    Parameters:\n",
    "    ====================\n",
    "    data: pandas dataframe containing columns of \"x\", \"y\" and covariates.\n",
    "    x_cols: list containing x to individually correlate to y\n",
    "    y_col: column header of y\n",
    "    covar: list of column headers of covariates\n",
    "    method: string, specify the method used for computing correlation coefficient, available methods are \"pearson\", \"spearman\", \"kandall\", \"shepherd\", etc. \n",
    "    Return:Stats pandas dataframe\n",
    "    ====================\n",
    "    \n",
    "    \"\"\"\n",
    "    df = data.copy()\n",
    "    columns_new = ['n', 'r', 'CI95%', 'p-val']\n",
    "    scores = pd.DataFrame(columns = columns_new)\n",
    "    for i in x_cols:\n",
    "        df = df.dropna(subset=[i, y_col])\n",
    "        result = pg.partial_corr(data=df, x=i, y=y_col, covar=covar, method=method)\n",
    "        result['protein'] = i\n",
    "        scores=scores.append(result, sort=False)\n",
    "    scores = scores.assign(new_column = lambda x: -np.log10(scores['p-val']))\n",
    "    scores = scores.rename({'new_column' : '-Log pvalue'}, axis=1)\n",
    "    \n",
    "    #FDR correction\n",
    "    np.random.seed(123)\n",
    "    reject, qvalue = multi.fdrcorrection(scores['p-val'], alpha=0.05, method='indep')\n",
    "    scores['qvalue'] = qvalue\n",
    "    scores['rejected'] = reject\n",
    "    scores = scores.set_index('protein')\n",
    "    return scores\n",
    "\n",
    "def return_partialcorr_sig(data, thresh=0.3):\n",
    "    \"\"\"\n",
    "    Takes output of partial_corr_pg, apply thresh and FDR corrected p-value as cutoff, return significant correlation matrix.\n",
    "    \"\"\"\n",
    "    df = data.copy()\n",
    "    sig = df.rejected & (abs(df.r)>=thresh)\n",
    "    df_sig = df[sig].sort_values(by='r', ascending=False)\n",
    "    return df_sig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experimental design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annotation_file = pd.read_csv('data/raw/Experiment annotation file_upgrade.csv', index_col = [0])\n",
    "annotation_file_plasma = annotation_file[annotation_file['Sample type'] == 'Plasma']\n",
    "annotation_file_liver = annotation_file[annotation_file['Sample type'] == 'Liver']\n",
    "annotation_file.groupby(['Sample type', 'Group2']).count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# ID mapping and external annotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IDmapping between sample ID, Participant ID and study groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IDmapping_samples = dict(zip(annotation_file_plasma['Participant ID'], annotation_file_plasma['Sample ID']))\n",
    "IDmapping_SampleID_to_Group2 = dict(zip(annotation_file_plasma['Sample ID'], annotation_file_plasma['Group2']))\n",
    "IDmapping_SampleID_to_ParticipantID=dict(zip(annotation_file['Sample ID'], annotation_file['Participant ID']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Protein functional/categorical annotations extracted from the Human Protein Atlas (HPA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tissue_specificity = pd.read_csv('data/annotations/HPA/tissue_specificity_rna_any.csv', index_col = [0])\n",
    "HPA_secreted=pd.read_csv('data/annotations/HPA/protein_class_Predicted_secreted.csv')\n",
    "secreted=HPA_secreted['Gene']\n",
    "secreted_dict={i:'secreted' for i in secreted}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ID mapping between lab test abbreviation and descriptions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labtest_id = pd.read_csv('data/raw/labtest_naming_key.csv')\n",
    "IDmapping_labtest_code_description = dict(zip(labtest_id['id_labtest'], labtest_id['description']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Clinical data processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Processing steps\n",
    "- Filter lab test items for 50% valid values across subjects in ald and hp cohort separately\n",
    "- Select overlapping lab test items in the ald and hp cohort\n",
    "- Merge ald and hp data into one dataframe\n",
    "- Biopsy-verified Kleiner score (a scoring system for staging liver fibrosis) in the high risk group spans between 0-4. In the at risk group, there was no appranet liver fibrosis assessed by FibroScan but no biopsy, this value was assigned to 'F0-1'. We convert 'F0-1'to 0.5\n",
    "- Convert Steatosis grade (<5%, 5-33%, 33-66%, >66%) to an ordinal score (0-3)\n",
    "- Convert gender, abstinent status and statin status to numerical values (male and yes for 1, otherwise 0)\n",
    "- Keep only numeric data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Clinical data for the healthy (hp) and disease (ald) cohort came in separate files, below pre-processing steps\n",
    "labtest_ald_raw = pd.read_csv('data/raw/GALA_ALD_all_lab_values.csv')\n",
    "labtest_hp_raw = pd.read_excel('data/raw/GALA_HP_proteomics_20180322_lab_values.xlsx')\n",
    "labtest_ald_raw = convert_to_numeric(labtest_ald_raw)\n",
    "labtest_hp_raw = convert_to_numeric(labtest_hp_raw)\n",
    "\n",
    "#Filter lab test for 50% valid values across subjects\n",
    "thresh = [data.shape[0] * 0.5 for data in [labtest_ald_raw, labtest_hp_raw]]\n",
    "labtest_ald = labtest_ald_raw.dropna(axis=1, thresh = thresh[0])\n",
    "labtest_hp = labtest_hp_raw.dropna(axis = 1, thresh = thresh[1])\n",
    "\n",
    "#Overlap of labtests between ald and hp\n",
    "labtest_overlap = list(set(labtest_hp.columns) & set(labtest_ald.columns))\n",
    "\n",
    "#Integrate lab tests\n",
    "labtest = labtest_ald.append(labtest_hp, ignore_index= True, sort = False)\n",
    "labtest['Sample ID'] = labtest['Participant ID'].map(IDmapping_samples)\n",
    "labtest = labtest.dropna(subset = ['Sample ID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = labtest[labtest.kleiner == 'F0-1'].index\n",
    "labtest.loc[mask,'kleiner'] = 0.5\n",
    "dict_nas_steatosis_ordinal = {value: counter for counter, value in enumerate(labtest.nas_steatosis.dropna().unique())}\n",
    "labtest['nas_steatosis_ordinal'] = labtest['nas_steatosis'].map(dict_nas_steatosis_ordinal)\n",
    "\n",
    "mask_hp = labtest[labtest['project'] == 'GALA-HP'].index\n",
    "for i in ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']:\n",
    "    labtest.loc[mask_hp, i] = -1\n",
    "\n",
    "labtest['gender_num'] = np.where(labtest['gender'].str.lower() == 'male', np.int64(1), np.int64(0))\n",
    "labtest['abstinent_num'] = np.where(labtest['abstinent'].str.lower() == 'yes', np.int64(1), np.int64(0))\n",
    "labtest['statin_num'] = np.where(labtest['medication_statin'].str.lower() == 'yes', np.int64(1), np.int64(0))\n",
    "\n",
    "labtest = labtest.fillna(np.nan)\n",
    "labtest = convert_to_numeric(labtest)\n",
    "labtest_numeric = labtest.set_index('Sample ID').select_dtypes(include = np.number).drop('id', axis = 1)\n",
    "\n",
    "labtest.to_csv('data/processed/labtest_integrated_all.csv')\n",
    "labtest_numeric.to_csv('data/processed/labtest_integrated_numeric.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Plasma proteomics data processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## dataset processing steps\n",
    "- Input data as Spectronaut v13 search engine output (details refer to Methods in the manuscript)\n",
    "- Filter at protein level for 60% data completeness across all runs\n",
    "- Filter at individual sample level for total quantified protein groups above 200\n",
    "- Log2 transform raw intensity values\n",
    "- Perform imputation on remaining missing values using a down-shifted Gaussian distribution across sample runs of the same protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report_plasma = pd.read_csv('data/raw/Protein_ALDupgrade_Report.csv', na_values='Filtered')\n",
    "report_plasma.rename({'PG.Genes': 'Gene names', \n",
    "                      'PG.ProteinAccessions': 'Protein IDs', \n",
    "                      'PG.ProteinDescriptions':'Protein names'}, inplace= True, axis=1)\n",
    "IDmapping_UniprotID_to_Genename=dict(zip(report_plasma['Protein IDs'], report_plasma['Gene names']))\n",
    "report_plasma=report_plasma.drop('Protein names', axis=1)\n",
    "experimental_columns = annotation_file_plasma['Sample ID']\n",
    "validation_columns = annotation_file_plasma[annotation_file_plasma['Groups']=='ALD-validation']['Sample ID'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_filenames_ids = dict(zip(annotation_file['File name'], annotation_file['Sample ID']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_raw = report_plasma.copy()\n",
    "data_plasma_raw.drop(data_plasma_raw.filter(regex='StrippedSequences').columns, axis=1, inplace = True)\n",
    "data_plasma_raw['Protein ID']=[i.split(';')[0] for i in data_plasma_raw['Protein IDs']]\n",
    "df = data_plasma_raw\n",
    "IDmapping_plasma_ProteinID_to_ProteinIDs = dict(zip(df['Protein ID'], df['Protein IDs']))\n",
    "data_plasma_raw = data_plasma_raw.rename(mapper = map_filenames_ids, axis=1).drop('Protein IDs', axis=1)\n",
    "data_plasma_raw.drop(validation_columns, axis=1, inplace=True)\n",
    "IDmapping_UniprotID_to_Genename = dict(zip(data_plasma_raw['Protein ID'], data_plasma_raw['Gene names']))\n",
    "data_plasma_raw.drop('Gene names', axis=1, inplace=True)\n",
    "data_plasma_raw.set_index('Protein ID', inplace=True)\n",
    "data_plasma_raw = data_plasma_raw.dropna(how='all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = data_plasma_raw.notna().sum(axis=1) > 603 * 0.6\n",
    "mask.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_COMPLETENESS = 0.6\n",
    "data_plasma_filtered = data_plasma_raw.dropna(axis=0, thresh = data_plasma_raw.shape[1] * DATA_COMPLETENESS)\n",
    "# data_plasma_filtered #hide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MIN_N_PROTEIN_GROUPS = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_filtered = data_plasma_filtered.dropna(axis=1, thresh = MIN_N_PROTEIN_GROUPS)\n",
    "# data_plasma_filtered #hide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_plasma_filtered.copy()\n",
    "data_plasma_filtered_log = log_transform_data(data.astype(np.float64), 2)\n",
    "data_plasma_filtered_log_imputed = imputation_normal_distribution(data_plasma_filtered_log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## data quality assessment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- QC samples are pooled plasma or liver homogenate from study samples, processed together with study samples as early as possible including sample preparation and LC-MS/MS analysis, thus indicate entire workflow variability across sample measurement period. \n",
    "- Take raw data as input dataset (before filtering/log2 transformation/imputation)\n",
    "- calculate coefficient of variation among all QC samples\n",
    "- Indicate number of proteins with a CV% < 30%, indicate median CV% among all proteins\n",
    "- Calculate pair-wise Pearson Correlation between QC samples "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data completeness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- before filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_plasma_raw.copy()\n",
    "data_completeness = data.count().sum()/(data.shape[0] * data.shape[1])\n",
    "print(\"Total: {}\".format(data.shape[0]))\n",
    "print(\"Median: {}\\nData completeness: {:.0%}\".format(data.count().median(), data_completeness))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- after filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_raw.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_plasma_filtered.copy()\n",
    "data_completeness = data.count().sum()/(data.shape[0] * data.shape[1])\n",
    "print(\"Total: {}\".format(data.shape[0]))\n",
    "print(\"Median: {}\\nData completeness: {:.0%}\".format(data.count().median(), data_completeness))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ploting figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qc_plasma = annotation_file_plasma[annotation_file_plasma['Group2'] == 'QC']['Sample ID']\n",
    "df_qc_plasma = data_plasma_raw.copy()\n",
    "df_qc_plasma = convert_to_numeric(df_qc_plasma)\n",
    "df_qc = df_qc_plasma[qc_plasma]\n",
    "df_qc = df_qc.dropna(thresh=2)\n",
    "func = lambda x: np.std(x) / np.mean(x)\n",
    "median = df_qc.median(axis = 1)\n",
    "var = df_qc.apply(func, axis = 1)\n",
    "df_qc = df_qc.assign(var = var, median = median)\n",
    "df_qc = df_qc.sort_values(by = 'median', ascending = False)\n",
    "df_qc = df_qc.assign(rank = np.arange(df_qc.shape[0]))\n",
    "df_qc['Gene name'] = df_qc.index.map(IDmapping_UniprotID_to_Genename)\n",
    "df_qc['median_log10'] = np.log10(df_qc['median'])\n",
    "df_qc_30 = df_qc[df_qc['var'] < 0.3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export source data for figS3b\n",
    "df_qc[['median_log10', 'var', 'rank', 'Gene name']].to_csv('data/processed/source_data_figS3b.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "fig, ax = plt.subplots(figsize = (4,4))\n",
    "plt.scatter(x = df_qc['median_log10'], y = df_qc['var'], facecolor = 'white', edgecolor = 'black')\n",
    "plt.scatter(x = df_qc_30['median_log10'], y = df_qc_30['var'], facecolor = 'white', edgecolor = 'darkblue')\n",
    "plt.ylabel('Coeff. of variation', fontsize = 15)\n",
    "plt.title('Total variability', fontsize = 15)\n",
    "plt.yticks(fontsize = 15, ticks = [0, 0.5, 1, 1.5, 2])\n",
    "plt.xticks(fontsize = 15)\n",
    "xlim = plt.gca().get_xlim()\n",
    "#plt.ylim(-0.1, 2)\n",
    "#plt.plot(xlim, [0.3, 0.3], color = 'black', lw = 1, ls = '--')\n",
    "plt.annotate(xy = [0.5, 0.9], xycoords = 'axes fraction', s = 'CV < 30%: {}'.format(df_qc_30.shape[0], fontsize = 12))\n",
    "plt.annotate(xy = [0.5, 0.85], xycoords = 'axes fraction', s = 'Median CV%: {}%'.format(df_qc['var'].median().round(2)*100), fontsize = 12)\n",
    "plt.xlabel('MS signal [Log10]', fontsize = 15)\n",
    "plt.savefig('figures/cv_plasma.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (4,4))\n",
    "plt.scatter(x = df_qc['rank'], y = df_qc['median_log10'], facecolor = 'white', edgecolor = 'black')\n",
    "plt.scatter(x = df_qc_30['rank'], y = df_qc_30['median_log10'], facecolor = 'white', edgecolor = 'darkblue')\n",
    "plt.ylabel('MS signal [Log10]', fontsize = 15)\n",
    "plt.title('Dynamic range', fontsize = 15)\n",
    "plt.xticks(fontsize = 15)\n",
    "plt.yticks(fontsize = 15)\n",
    "xlim = plt.gca().get_xlim()\n",
    "plt.ylim(2.5, 9.2)\n",
    "plt.xlabel('Abundance rank', fontsize = 15)\n",
    "plt.savefig('figures/dynamic range_plasma.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_filtered_log[qc_plasma].corr().to_csv('data/processed/source_data_figS3a.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (6,5))\n",
    "sns.set_style('ticks')\n",
    "sns.heatmap(data_plasma_filtered_log[qc_plasma].corr(), \n",
    "            vmin = 0.5, vmax = 1, lw =2, annot = False, cbar = True,\n",
    "           cmap = 'Blues')\n",
    "sns.set(font_scale=0.9)\n",
    "plt.title('Pearson correlation', fontsize = 15)\n",
    "plt.yticks(fontsize = 12)\n",
    "plt.xticks(ticks = None, fontsize = 12)\n",
    "plt.savefig('figures/qc_corr_plasma.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Liver proteomics data processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset processing\n",
    "- Input data as Spectronaut v13 search engine output (details refer to Methods in the manuscript)\n",
    "- Filter at protein level for 60% data completeness across all runs\n",
    "- Filter at individual sample level for total quantified protein groups above 4000\n",
    "- Log2 transform raw intensity values\n",
    "- Perform imputation on remaining missing values using a down-shifted Gaussian distribution across sample runs of the same protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report_liver = pd.read_csv('data/raw/Protein_20200221_121354_20200218_ALD_LiverTissue_PlateS1_Atlaslib_Report_PGonly.csv',\n",
    "                          na_values='Filtered')\n",
    "data_liver_raw = report_liver.copy()\n",
    "data_liver_raw = data_liver_raw.rename({'PG.ProteinAccessions': 'Protein IDs', 'PG.Genes':'Gene names'}, axis=1)\n",
    "data_liver_raw['Protein ID']=[i.split(';')[0] for i in data_liver_raw['Protein IDs']]\n",
    "data_liver_raw['Gene name'] = [str(i).split(';')[0] for i in data_liver_raw['Gene names']]\n",
    "df=data_liver_raw\n",
    "IDmapping_liver_ProteinID_to_ProteinIDs = dict(zip(df['Protein ID'], df['Protein IDs']))\n",
    "IDmapping_Liver_UniprotID_to_Genename = dict(zip(df['Protein ID'], df['Gene name']))\n",
    "to_drop = ['PG.Qvalue', 'PG.MolecularWeight', 'PG.ProteinDescriptions', 'Protein IDs', 'Gene names', 'Gene name']\n",
    "\n",
    "data_liver_raw = data_liver_raw.rename(mapper=map_filenames_ids, axis=1).drop(to_drop, axis=1).set_index('Protein ID')\n",
    "data_liver_raw = data_liver_raw.dropna(how='all')\n",
    "\n",
    "data_liver_filtered = data_liver_raw.dropna(thresh = data_liver_raw.shape[1] * 0.6)\n",
    "columns_count = [i for i in data_liver_filtered.columns if data_liver_filtered[i].count() > 4000]\n",
    "data_liver_filtered = data_liver_filtered[columns_count]\n",
    "data_liver_filtered = convert_to_numeric(data_liver_filtered)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = data_liver_filtered.copy()\n",
    "data_liver_filtered_log = log_transform_data(df.astype(np.float64), 2)\n",
    "data_liver_filtered_log_imputed = imputation_normal_distribution(data_liver_filtered_log)\n",
    "data_liver_filtered_log_imputed = data_liver_filtered_log_imputed.rename_axis('Sample ID', axis =1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_RESULTS = 'data/processed/proteomics datasets.xlsx'\n",
    "\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    data_liver_raw.to_excel(writer, sheet_name='data_liver_raw')\n",
    "    data_liver_filtered.to_excel(writer, sheet_name='data_liver_filtered')\n",
    "    data_liver_filtered_log_imputed.to_excel(writer, sheet_name='data_liver_filtered_log2_imputed')\n",
    "    data_plasma_raw.to_excel(writer, sheet_name='data_plasma_raw')\n",
    "    data_plasma_filtered.to_excel(writer, sheet_name='data_plasma_filtered')\n",
    "    data_plasma_filtered_log_imputed.to_excel(writer, sheet_name='data_plasma_filtered_log2_imputed') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data quality assessment\n",
    "- QC samples are pooled plasma or liver homogenate from study samples, processed together with study samples as early as possible including sample preparation and LC-MS/MS analysis, thus indicate entire workflow variability across sample measurement period. \n",
    "- Take raw data as input dataset (before filtering/log2 transformation/imputation)\n",
    "- calculate coefficient of variation among all QC samples\n",
    "- Indicate number of proteins with a CV% < 30%, indicate median CV% among all proteins\n",
    "- Calculate pair-wise Pearson Correlation between QC samples "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qc_liver = annotation_file_liver[annotation_file_liver['Groups'] == 'QC_liver']['Sample ID'].to_list()\n",
    "df_qc_liver = data_liver_raw.copy()\n",
    "df_qc_liver = convert_to_numeric(df_qc_liver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "df_qc = df_qc_liver[qc_liver].dropna(thresh = 2)\n",
    "func = lambda x: np.std(x) / np.mean(x)\n",
    "median = df_qc.median(axis = 1)\n",
    "var = df_qc.apply(func, axis = 1)\n",
    "df_qc = df_qc.assign(var = var, median = median)\n",
    "df_qc = df_qc.sort_values(by = 'median', ascending = False)\n",
    "df_qc = df_qc.assign(rank = np.arange(df_qc.shape[0]))\n",
    "df_qc['Gene name'] = df_qc.index.map(IDmapping_Liver_UniprotID_to_Genename)\n",
    "df_qc['median_log10'] = np.log10(df_qc['median'])\n",
    "df_qc_30 = df_qc[df_qc['var'] < 0.3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export source data for figS3b\n",
    "df_qc[['median_log10', 'var', 'rank', 'Gene name']].to_csv('data/processed/source_data_figS3e.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (4,4))\n",
    "plt.scatter(x = df_qc['median_log10'], y = df_qc['var'], facecolor = 'white', edgecolor = 'black')\n",
    "plt.scatter(x = df_qc_30['median_log10'], y = df_qc_30['var'], facecolor = 'white', edgecolor = 'darkblue')\n",
    "plt.ylabel('Coeff. of variation', fontsize = 15)\n",
    "plt.title('Total variability', fontsize = 15)\n",
    "plt.yticks(fontsize = 15, ticks = [0, 0.5, 1, 1.5, 2])\n",
    "plt.xticks(fontsize = 15)\n",
    "xlim = plt.gca().get_xlim()\n",
    "plt.ylim(-0.1, 2)\n",
    "#plt.plot(xlim, [0.3, 0.3], color = 'black', lw = 1, ls = '--')\n",
    "plt.annotate(xy = [0.5, 0.9], xycoords = 'axes fraction', s = 'CV < 30%: {}'.format(df_qc_30.shape[0], fontsize = 12))\n",
    "plt.annotate(xy = [0.5, 0.85], xycoords = 'axes fraction', s = 'Median CV%: {}%'.format(df_qc['var'].median().round(2)*100), fontsize = 12)\n",
    "plt.xlabel('MS signal [Log10]', fontsize = 15)\n",
    "plt.savefig('figures/cv_liver.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "fig, ax = plt.subplots(figsize = (4,4))\n",
    "plt.scatter(x = df_qc['rank'], y = df_qc['median_log10'], facecolor = 'white', edgecolor = 'black')\n",
    "plt.scatter(x = df_qc_30['rank'], y = df_qc_30['median_log10'], facecolor = 'white', edgecolor = 'darkblue')\n",
    "plt.ylabel('MS signal [Log10]', fontsize = 15)\n",
    "plt.title('Dynamic range', fontsize = 15)\n",
    "plt.xticks(fontsize = 15)\n",
    "plt.yticks(fontsize = 15)\n",
    "xlim = plt.gca().get_xlim()\n",
    "plt.ylim(2.5, 9.2)\n",
    "plt.xlabel('Abundance rank', fontsize = 15)\n",
    "plt.savefig('figures/dynamicRange_liver.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_liver_filtered_log[qc_liver].corr().to_csv('data/processed/source_data_figS3d.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize = (6,5))\n",
    "sns.set_style('ticks')\n",
    "sns.heatmap(data_liver_filtered_log[qc_liver].corr(), \n",
    "            vmin = 0.5, vmax = 1, lw =2, annot = False, cbar = True,\n",
    "           cmap='Blues')\n",
    "sns.set(font_scale=0.9)\n",
    "plt.title('Pearson correlation', fontsize = 15)\n",
    "plt.yticks(fontsize = 12)\n",
    "plt.xticks(ticks = None, fontsize = 12)\n",
    "plt.savefig('figures/qc_corr_liver.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connecting clinical data to liver proteome data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cli_liver = labtest_numeric.copy()\n",
    "df_cli_liver['Participant ID'] = df_cli_liver.index.map(dict(zip(annotation_file_plasma['Sample ID'], annotation_file_plasma['Participant ID'])))\n",
    "df_cli_liver['Sample ID liver'] = df_cli_liver['Participant ID'].map(dict(zip(annotation_file_liver['Participant ID'], annotation_file_liver['Sample ID'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cli_liver = df_cli_liver.dropna(subset = ['Sample ID liver']).set_index('Sample ID liver')\n",
    "IDmapping_liverSampleID_to_ParticipantID = dict(zip(df_cli_liver.index, df_cli_liver['Participant ID']))\n",
    "para_tokeep = ['age', 'kleiner', 'nas', 'nas_inflam', 'nas_steatosis_ordinal', 'gender_num', 'bmi', 'abstinent_num']\n",
    "dfliver_perseus = data_liver_filtered_log_imputed.copy().T\n",
    "dfliver_perseus = dfliver_perseus.join(df_cli_liver[para_tokeep]).T\n",
    "dfliver_perseus['Gene names'] = dfliver_perseus.index.map(IDmapping_Liver_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# ANCOVA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Preparing dataframes (wide- to long format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wide_to_long(df_prot, df_cli, para_to_combine, sample_type='liver'):\n",
    "    data=df_prot.copy().rename_axis('Sample ID', axis=1)\n",
    "    data_long=pd.melt(data, value_vars=list(data.columns), value_name='Intensity')\n",
    "    data_long['Protein ID']=np.tile(data.index, data.shape[1])\n",
    "    data_long=data_long.assign(type=lambda x: sample_type).set_index('Sample ID')\n",
    "    data_long_new=data_long.join(df_cli[para_to_combine]).reset_index().set_index('Protein ID')\n",
    "    return data_long_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "para_to_combine=['age', 'kleiner', 'nas', 'nas_inflam', 'nas_steatosis_ordinal', 'gender_num', 'bmi', 'abstinent_num']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Long format datasets with unimptuted proteomics data\n",
    "- Liver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prot=data_liver_filtered_log.copy()\n",
    "df_cli=df_cli_liver.copy()\n",
    "data_liver_filtered_log_long=wide_to_long(df_prot=df_prot, df_cli=df_cli, para_to_combine=para_to_combine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prot=data_plasma_filtered_log.copy()\n",
    "df_cli=labtest_numeric.copy()\n",
    "data_plasma_filtered_log_long=wide_to_long(df_prot=df_prot, df_cli=df_cli, para_to_combine=para_to_combine, sample_type='plasma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Long format datasets with imputed proteomics data\n",
    "- Liver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prot=data_liver_filtered_log_imputed.copy()\n",
    "df_cli=df_cli_liver.copy()\n",
    "data_liver_filtered_log_imputed_long=wide_to_long(df_prot=df_prot, df_cli=df_cli, para_to_combine=para_to_combine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_liver_rename = data_liver_raw.copy().rename(mapper=dict(zip(IDmapping_liverSampleID_to_ParticipantID.values(),IDmapping_liverSampleID_to_ParticipantID.keys())), axis=1)\n",
    "data_liver_long = wide_to_long(df_prot=data_liver_rename, df_cli=df_cli_liver, para_to_combine=para_to_combine, sample_type='liver')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prot=data_plasma_filtered_log_imputed.copy()\n",
    "df_cli=labtest_numeric.copy()\n",
    "data_plasma_filtered_log_imputed_long=wide_to_long(df_prot=df_prot, df_cli=df_cli, para_to_combine=para_to_combine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Executing ANCOVA in the plasma proteome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Covariates included in ANCOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp = annotation_file_plasma[annotation_file_plasma['Group2'] == 'HP']['Sample ID']\n",
    "data = data_plasma_filtered_log_imputed_long.copy()\n",
    "data = data[np.logical_not(data['Sample ID'].isin(hp))]\n",
    "data_fibrosis = data[data['kleiner']!= 0.5]\n",
    "covar = ['age', 'bmi', 'gender_num', 'abstinent_num', 'nas_steatosis_ordinal']\n",
    "covar_steatosis = ['age', 'bmi', 'gender_num', 'abstinent_num', 'kleiner', 'nas_inflam']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RECALCULATE_ANCOVA = True\n",
    "if not RECALCULATE_ANCOVA:\n",
    "    ancova_results_file = os.path.join('data/processed', 'ANCOVA_results.xlsx')\n",
    "    result_ancova_kleiner_p = pd.read_excel(ancova_results_file, sheet_name='ANCOVA_plasma_Kleiner')\n",
    "    result_ancova_inflam_p = pd.read_excel(ancova_results_file, sheet_name='ANCOVA_plasma_inflammation')\n",
    "    result_ancova_steatosis_p = pd.read_excel(ancova_results_file, sheet_name='ANCOVA_plasma_steatosis')\n",
    "else: \n",
    "    result_ancova_kleiner_p = ancova_pg(data = data_fibrosis, dv = 'Intensity', between = 'kleiner', covar = covar)\n",
    "    result_ancova_inflam_p = ancova_pg(data = data, dv = 'Intensity', between = 'nas_inflam', covar = covar)\n",
    "    result_ancova_steatosis_p = ancova_pg(data = data, dv = 'Intensity', between = 'nas_steatosis_ordinal', covar = covar_steatosis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_sig_ancova(ancova_result, source, mapper=IDmapping_UniprotID_to_Genename):\n",
    "    df = ancova_result.copy()\n",
    "    sig_df = df[(df.rejected) & (df['Source']==source)]\n",
    "    sig_df.loc[:, 'Gene name'] = sig_df['protein'].map(mapper)\n",
    "    return (sig_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract significant proteins from ANCOVA results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_ancova_kleiner_p = extract_sig_ancova(result_ancova_kleiner_p, 'kleiner')\n",
    "sig_ancova_inflam_p = extract_sig_ancova(result_ancova_inflam_p, 'nas_inflam')\n",
    "sig_ancova_steatosis_p = extract_sig_ancova(result_ancova_steatosis_p, 'nas_steatosis_ordinal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Venn diagram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib_venn import venn3\n",
    "fig, ax = plt.subplots(figsize = (4, 4))\n",
    "sig_list_p = [sig_ancova_kleiner_p.protein, sig_ancova_inflam_p.protein, sig_ancova_steatosis_p.protein]\n",
    "v = venn3(subsets = [set(i) for i in sig_list_p], set_labels = ['fibrosis', 'inflammation', 'steatosis'], alpha = 0.5,\n",
    "         set_colors=('steelblue', 'g', 'lightyellow'))\n",
    "plt.savefig('figures/ancova_plasma_venn.pdf', dpi = 80, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract significant proteins ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fib_p = sig_ancova_kleiner_p['protein'].to_list()\n",
    "inf_p = sig_ancova_inflam_p['protein'].to_list()\n",
    "stea_p = sig_ancova_steatosis_p['protein'].to_list()\n",
    "ancova_all_p = list(set(fib_p + inf_p + stea_p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ANCOVA in plasma poteome on the subset of n=79 where liver proteomes were measured"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_79_participantID = annotation_file[(annotation_file['Sample type']=='Liver') & (annotation_file['Groups']!='QC_liver')]['Participant ID']\n",
    "subset_79_plasma_sampleID = annotation_file[annotation_file['Sample type']=='Plasma'].set_index('Participant ID').loc[subset_79_participantID]['Sample ID']\n",
    "data_plasma_79 = data[data['Sample ID'].isin(subset_79_plasma_sampleID)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ancova_79_folder = 'data/processed/ancova_79_plasma'\n",
    "if not RECALCULATE_ANCOVA:\n",
    "    ancova_k_p_79 = pd.read_pickle(os.path.join(ancova_79_folder, 'ancova_79_fibrosis.pkl'))\n",
    "    ancova_i_p_79 = pd.read_pickle(os.path.join(ancova_79_folder, 'ancova_79_inflam.pkl'))\n",
    "    ancova_s_p_79 = pd.read_pickle(os.path.join(ancova_79_folder, 'ancova_79_steatosis.pkl'))\n",
    "else:\n",
    "    ancova_k_p_79 = ancova_pg(data = data_plasma_79, dv = 'Intensity', between = 'kleiner', covar = covar)\n",
    "    ancova_i_p_79 = ancova_pg(data = data_plasma_79, dv = 'Intensity', between = 'nas_inflam', covar = covar)\n",
    "    ancova_s_p_79 = ancova_pg(data = data_plasma_79, dv = 'Intensity', between = 'nas_steatosis_ordinal', covar = covar_steatosis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ancova_79_f = extract_sig_ancova(ancova_k_p_79, 'kleiner')\n",
    "ancova_79_i = extract_sig_ancova(ancova_i_p_79, 'nas_inflam')\n",
    "ancova_79_s = extract_sig_ancova(ancova_s_p_79, 'nas_steatosis_ordinal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancova_79_all = set(ancova_79_f['protein'].tolist() + ancova_79_i['protein'].tolist() + ancova_79_s['protein'].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(ancova_79_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Preparing supplementary tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_plasma_79.copy()\n",
    "df_mean = data.groupby([data.index, 'kleiner']).mean()['Intensity'].unstack()\n",
    "df_mean.columns=['F' + str(i) for i in range(5)]\n",
    "for i in [2, 3, 4]:\n",
    "    col_name = 'F{}/F0'.format(i)\n",
    "    df_mean[col_name] = 2** (df_mean['F{}'.format(i)]-df_mean['F0'])\n",
    "df_mean['Gene name']=df_mean.index.map(IDmapping_UniprotID_to_Genename)\n",
    "col_names = ['ANCOVA significant '+i for i in ['fibrosis', 'inflammation', 'steatosis']]\n",
    "annotations = [i['protein'].tolist() for i in [ancova_79_f, ancova_79_i, ancova_79_s]]\n",
    "for col_name, annotation in zip(col_names, annotations):\n",
    "    df_mean[col_name]= np.where(df_mean.index.isin(annotation), '+', '')\n",
    "df_mean['ANCOVA significant combined'] = np.where(df_mean.index.isin(ancova_79_all), '+', '')\n",
    "df_mean['Protein IDs']=df_mean.index.map(IDmapping_plasma_ProteinID_to_ProteinIDs)\n",
    "data_mean_plasma_subset = df_mean.copy()\n",
    "df_mean.to_csv('data/processed/ancova_79_plasma/ANCOVA_plasma_79subset_sig.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Executing ANCOVA in the liver proteome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_liver_filtered_log_imputed_long.copy()\n",
    "if not RECALCULATE_ANCOVA:\n",
    "    ancova_results_file = os.path.join('data/processed', 'ANCOVA_results.xlsx')\n",
    "    result_ancova_kleiner_l = pd.read_excel(ancova_results_file, sheet_name='ANCOVA_liver_Kleiner')\n",
    "    result_ancova_inflam_l = pd.read_excel(ancova_results_file, sheet_name='ANCOVA_liver_inflammation')\n",
    "    result_ancova_steatosis_l = pd.read_excel(ancova_results_file, sheet_name='ANCOVA_liver_steatosis')\n",
    "else: \n",
    "    result_ancova_kleiner_l = ancova_pg(data = data, dv = 'Intensity', between = 'kleiner', covar = covar)\n",
    "    result_ancova_inflam_l = ancova_pg(data = data, dv = 'Intensity', between = 'nas_inflam', covar = covar)\n",
    "    result_ancova_steatosis_l = ancova_pg(data = data, dv = 'Intensity', between = 'nas_steatosis_ordinal', covar = covar_steatosis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract significant proteins from ANCOVA results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapper_liver = IDmapping_Liver_UniprotID_to_Genename\n",
    "sig_ancova_kleiner_l = extract_sig_ancova(result_ancova_kleiner_l, 'kleiner', mapper=mapper_liver)\n",
    "sig_ancova_inflam_l = extract_sig_ancova(result_ancova_inflam_l, 'nas_inflam', mapper=mapper_liver)\n",
    "sig_ancova_steatosis_l = extract_sig_ancova(result_ancova_steatosis_l, 'nas_steatosis_ordinal', mapper=mapper_liver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting significant protein ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fib_l = sig_ancova_kleiner_l['protein'].to_list()\n",
    "inf_l = sig_ancova_inflam_l['protein'].to_list()\n",
    "stea_l = sig_ancova_steatosis_l['protein'].to_list()\n",
    "ancova_all_l = list(set(fib_l + inf_l + stea_l))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Venn diagram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib_venn import venn3\n",
    "fig, ax = plt.subplots(figsize = (4, 4))\n",
    "sig_list_l = [sig_ancova_kleiner_l.protein, sig_ancova_inflam_l.protein, sig_ancova_steatosis_l.protein]\n",
    "v = venn3(subsets = [set(i) for i in sig_list_l], set_labels = ['fibrosis', 'inflammation', 'steatosis'], alpha = 0.5,\n",
    "         set_colors=('steelblue', 'g', 'lightyellow'))\n",
    "plt.savefig('figures/ancova_liver_venn.pdf', dpi = 80, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_RESULTS = 'data/processed/ANCOVA_results.xlsx'\n",
    "\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    result_ancova_kleiner_p.to_excel(writer, sheet_name='ANCOVA_plasma_Kleiner')\n",
    "    result_ancova_inflam_p.to_excel(writer, sheet_name='ANCOVA_plasma_inflammation')\n",
    "    result_ancova_steatosis_p.to_excel(writer, sheet_name='ANCOVA_plasma_steatosis')\n",
    "   \n",
    "    result_ancova_kleiner_l.to_excel(writer, sheet_name='ANCOVA_liver_Kleiner')\n",
    "    result_ancova_inflam_l.to_excel(writer, sheet_name='ANCOVA_liver_inflammation')\n",
    "    result_ancova_steatosis_l.to_excel(writer, sheet_name='ANCOVA_liver_steatosis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing results as supplementary tables "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Liver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=data_liver_filtered_log_long.copy()\n",
    "df_mean = data.groupby([data.index, 'kleiner']).mean()['Intensity'].unstack()\n",
    "df_mean.columns=['F'+ str(i) for i in range(5)]\n",
    "for i in range(2,5):\n",
    "    df_mean['F{}/F0'.format(i)] = 2**(df_mean['F{}'.format(i)] - df_mean['F0'])\n",
    "\n",
    "df_mean['Gene name'] = df_mean.index.map(IDmapping_Liver_UniprotID_to_Genename)\n",
    "items = ['ANCOVA significant liver ' + item for item in ['combined', 'fibrosis', 'inflammation', 'steatosis']]\n",
    "df_items = [ancova_all_l, fib_l, inf_l, stea_l]\n",
    "for i, j in zip(items, df_items):\n",
    "    df_mean[i] = np.where(df_mean.index.isin(j), '+', '')\n",
    "    \n",
    "df_mean = df_mean.sort_values(by = 'F4/F0', ascending = False)\n",
    "df_mean['Protein IDs']=df_mean.index.map(IDmapping_liver_ProteinID_to_ProteinIDs)\n",
    "data_liver_mean_kleiner = df_mean.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_liver_mean_kleiner['ANCOVA significant liver combined'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_plasma_filtered_log_long.copy()\n",
    "df_mean = data.groupby([data.index, 'kleiner']).mean()['Intensity'].unstack()\n",
    "df_mean = df_mean.drop(0.5, axis = 1)\n",
    "df_mean.columns=['HP'] + ['F'+ str(i) for i in range(5)]\n",
    "for i in range(2,5):\n",
    "    df_mean['F{}/F0'.format(i)] = 2**(df_mean['F{}'.format(i)] - df_mean['F0'])\n",
    "    \n",
    "df_mean['Gene name'] = df_mean.index.map(IDmapping_UniprotID_to_Genename)\n",
    "items = ['ANCOVA significant plasma ' + item for item in ['combined', 'fibrosis', 'inflammation', 'steatosis']]\n",
    "df_items = [ancova_all_p, fib_p, inf_p, stea_p]\n",
    "for i, j in zip(items, df_items):\n",
    "    df_mean[i] = np.where(df_mean.index.isin(j), '+', '')\n",
    "    \n",
    "df_mean = df_mean.sort_values(by = 'F4/F0', ascending = False)\n",
    "df_mean['Protein IDs']=df_mean.index.map(IDmapping_plasma_ProteinID_to_ProteinIDs)\n",
    "data_plasma_mean_kleiner = df_mean.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_mean_kleiner['ANCOVA significant plasma combined'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save results (significant proteins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_RESULTS = 'data/processed/ANCOVA_SIG_GroupMeans.xlsx'\n",
    "\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    data_liver_mean_kleiner.to_excel(writer, sheet_name='ANCOVA_liver_mean')\n",
    "    data_plasma_mean_kleiner.to_excel(writer, sheet_name='ANCOVA_plasma_mean')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare for supplementary figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1=data_plasma_mean_kleiner.copy()\n",
    "df1=df1[df1['ANCOVA significant plasma combined']=='+']\n",
    "df1=df1[['F4/F0']]\n",
    "df2=data_mean_plasma_subset.copy()\n",
    "df2=df2[df2['ANCOVA significant combined']=='+']\n",
    "df3 = df1[['F4/F0']].join(df2[['F4/F0']].rename({'F4/F0':'F4/F0 subset'}, axis=1))\n",
    "df3.to_csv('data/processed/source_data_figS5f.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax=plt.subplots(figsize=(20,6))\n",
    "plt.bar(x=np.arange(df3.shape[0]), height=df3['F4/F0'], color='lightblue', )\n",
    "plt.bar(x=np.arange(df3.shape[0]), height=df3['F4/F0 subset'], color='darkred')\n",
    "plt.ylabel('F4/F0', fontsize=16)\n",
    "plt.xlabel('Protein rank based on F4/F0', fontsize=16)\n",
    "plt.savefig('figures/samplesize_effectsize.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### mapping co-dysregulated proteins to dysregulated proteins in plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_co = pd.read_csv('data/processed/co-dysregulated.csv')\n",
    "co=df_co['plasma'].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in co:\n",
    "    df3.loc[i, 'co-dysregulated']=df3.loc[i, 'F4/F0']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax=plt.subplots(figsize=(20,6))\n",
    "plt.bar(x=np.arange(df3.shape[0]), height=df3['F4/F0'], color='lightblue', )\n",
    "plt.bar(x=np.arange(df3.shape[0]), height=df3['co-dysregulated'], color='darkred')\n",
    "plt.ylabel('F4/F0', fontsize=16)\n",
    "plt.xlabel('Protein rank based on F4/F0', fontsize=16)\n",
    "plt.savefig('figures/co-dysregulated_to_all_plasma.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df3[df3['F4/F0']>1]['co-dysregulated'].count()/df3[df3['F4/F0']>1].shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.statistical_testing import pg_ttest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=data_liver_filtered_log_imputed_long.copy()\n",
    "ttest_liver_f4vsf0=pg_ttest(data=data, group_col='kleiner', group1=0, group2=4, value_col='Intensity')\n",
    "data=data_plasma_filtered_log_imputed_long.copy()\n",
    "ttest_plasma_f4vsf0=pg_ttest(data=data, group_col='kleiner', group1=0, group2=4, value_col='Intensity')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_liver_f4f0 = ttest_liver_f4vsf0[ttest_liver_f4vsf0['rejected']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax=plt.subplots(figsize=(6,6))\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.scatterplot(x='difference', y='-Log pvalue', data=ttest_liver_f4vsf0, hue='rejected')\n",
    "plt.xlim(-4.5, 4.5)\n",
    "plt.savefig('figures/volcano_f4f0.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax=plt.subplots(figsize=(6,6))\n",
    "ttest_liver_f4vsf0['difference'].hist(bins=40)\n",
    "sig_liver_f4f0['difference'].hist(bins=20)\n",
    "plt.savefig('figures/histgram.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Integration of liver- and plasma proteome data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "## Prepare dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma = data_plasma_raw.copy()\n",
    "data_liver = data_liver_raw.copy()\n",
    "print('plasma dataset: {}\\nliver dataset:{}'.format(data_plasma.shape, data_liver.shape))\n",
    "\n",
    "data_plasma = data_plasma.rename(mapper = dict(zip(annotation_file_plasma['Sample ID'], annotation_file_plasma['Participant ID'])), axis = 1)\n",
    "data_liver = data_liver.rename(mapper = dict(zip(annotation_file_liver['Sample ID'], annotation_file_liver['Participant ID'])), axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlap_index = pd.read_csv('data/processed/overlap_420_upgrade.csv')\n",
    "overlap_liver = overlap_index['liver']\n",
    "overlap_plasma = overlap_index['plasma']\n",
    "dict_plasma_IDswap = dict(zip(overlap_plasma.tolist(), overlap_liver.tolist()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlap_subjects = list(set(data_plasma.columns) & set(data_liver.columns))\n",
    "print('overlap subject ID number: {}'.format(len(overlap_subjects)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_overlap = data_plasma.loc[overlap_plasma, overlap_subjects]\n",
    "data_plasma_overlap = data_plasma_overlap.rename(mapper = dict_plasma_IDswap, axis=0)\n",
    "data_liver_overlap = data_liver.loc[overlap_liver, overlap_subjects]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_patiID_kleiner = dict(zip(labtest['Participant ID'], labtest['kleiner']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Protein level in liver-plasma space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = data_plasma_raw.copy().loc[overlap_plasma].T.median()\n",
    "df2 = data_liver_raw.copy().loc[overlap_liver].T.median()\n",
    "df1.index = df1.index.map(dict_plasma_IDswap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "liver_plasma = pd.concat([df1, df2], axis=1, )\n",
    "liver_plasma.columns=['plasma_median', 'liver_median']\n",
    "liver_plasma.to_csv('data/processed/liver_plasma_median_420.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting data to long format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "long_plasma = pd.melt(data_plasma_overlap, value_vars = list(data_plasma_overlap.columns), value_name = 'Intensity')\n",
    "long_plasma['index'] = np.tile(data_plasma_overlap.index, data_plasma_overlap.shape[1])\n",
    "long_plasma = long_plasma.assign(type = lambda x : 'plasma')\n",
    "\n",
    "long_liver = pd.melt(data_liver_overlap, value_vars = list(data_liver_overlap.columns), value_name = 'Intensity')\n",
    "long_liver['index'] = np.tile(data_liver_overlap.index, data_liver_overlap.shape[1])\n",
    "long_liver = long_liver.assign(type = lambda x : 'liver')\n",
    "\n",
    "long_combined = pd.merge(left= long_plasma, right = long_liver, left_on = ['variable', 'index'], right_on = ['variable', 'index'])\n",
    "long_combined['kleiner'] = long_combined['variable'].map(dict_patiID_kleiner)\n",
    "long_combined = long_combined.set_index('index').rename_axis('Protein ID', axis= 0)\n",
    "long_combined = long_combined.rename({'variable':'Participant ID'})\n",
    "long_combined = convert_to_numeric(long_combined)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exporting data for Dash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_overlap_median = pd.DataFrame(data_plasma_overlap.T.median(), columns=['plasma'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PLoverlap_median = pd.DataFrame(data_liver_overlap.T.median(), columns=['liver']).merge(pd.DataFrame(data_plasma_overlap.T.median(), columns=['plasma']), on='Protein ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PLoverlap_median = pd.DataFrame(data_liver_overlap.T.median(), columns=['liver']).merge(pd.DataFrame(data_plasma_overlap.T.median(), columns=['plasma']), on='Protein ID')\n",
    "PLoverlap_median = np.log10(PLoverlap_median)\n",
    "long_combined_export=long_combined.copy()\n",
    "for i in [PLoverlap_median, long_combined_export]:\n",
    "    i['Gene name'] = i.index.map(IDmapping_UniprotID_to_Genename)\n",
    "    i['Gene name'] = np.where(i['Gene name'].isnull(), i.index, i['Gene name'])\n",
    "    i['Genename_ProteinID']=i['Gene name']+'__'+i.index\n",
    "\n",
    "long_combined_export.rename({'Intensity_x':'plasma', 'Intensity_y':'liver'}, axis=1, inplace=True)\n",
    "\n",
    "for col in ['liver', 'plasma']:\n",
    "    long_combined_export[col]=np.log10(long_combined_export[col])\n",
    "    \n",
    "long_combined_export.dropna(subset=['plasma', 'liver']).drop(['variable'], axis=1).to_csv('../ALD-App/datasets/PLoverlap_data.csv')\n",
    "PLoverlap_median.dropna(subset=['liver', 'plasma']).to_csv('../ALD-App/datasets/PLoverlap_median.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "long_combined_export"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_corr = pd.merge(long_liver, long_plasma, left_on = ['variable', 'index'], right_on = ['variable', 'index'])\n",
    "data_corr.columns = ['Participant ID', 'Intensity_liver', 'Protein ID', 'type_liver', 'Intensity_plasma', 'type_plasma']\n",
    "data_corr = pd.merge(data_corr, df_cli_liver, left_on = 'Participant ID', right_on = 'Participant ID' )\n",
    "data_corr = convert_to_numeric(data_corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Compute pairwise correlation between liver and plasma protein levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plasma_liver_corr = pairwise_correlation(data = long_combined.drop(['kleiner'], axis = 1), columns = ['Intensity_x'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add HPA annotations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plcorr_tissue = pd.merge(left = plasma_liver_corr, right = tissue_specificity[['RNA tissue category', 'RNA TS TPM']], \n",
    "                         left_on = 'Gene name', right_on = 'Gene', how = 'left')\n",
    "plcorr_tissue.index = plasma_liver_corr.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_false = plcorr_tissue[plcorr_tissue['rejected'] == False].index\n",
    "corr_true = plcorr_tissue[(plcorr_tissue.rejected) & (abs(plcorr_tissue['r']) >= 0.3)].index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting significant correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_corr = plcorr_tissue.loc[corr_true].sort_values(by = 'r', ascending = False)\n",
    "sig_corr_top= sig_corr[:20].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_corr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig_corr.to_csv('data/processed/plasma_liver_corr_sig.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_plot_corr(protein, df):\n",
    "    df_protein = df.loc[protein].dropna(subset=['Intensity_liver', 'Intensity_plasma'])\n",
    "    #export source data\n",
    "    df_protein_source = df_protein[['Intensity_liver', 'Intensity_plasma', 'kleiner', 'nas_inflam']]\n",
    "    df_protein_source.to_csv('data/processed/source_data_fig3d_{}.csv'.format(protein))\n",
    "    r = str(plasma_liver_corr.loc[protein]['r'])\n",
    "    color_scheme='kleiner'\n",
    "    genename = IDmapping_Liver_UniprotID_to_Genename[protein]\n",
    "    fig, ax = plt.subplots(figsize= (4,4))\n",
    "    x = df_protein['Intensity_plasma']\n",
    "    y = df_protein['Intensity_liver']\n",
    "    x = np.log2(x)\n",
    "    y = np.log2(y)\n",
    "    plt.scatter(x, y, alpha = 1, s= 140, c= df_protein[color_scheme], cmap = 'coolwarm', edgecolor = 'black')\n",
    "    plt.xlabel('Intensity in plasma [Log2]', fontsize = 15)\n",
    "    plt.xticks(fontsize=15)\n",
    "    plt.yticks(fontsize=15)\n",
    "    plt.ylabel('Intensity in liver [Log2]', fontsize = 15)\n",
    "    plt.title(genename, fontsize = 30)\n",
    "    plt.text(x=0.05, y=0.9, transform=ax.transAxes, s= 'r = {:.2f}'.format(float(r)), fontsize = 25, color='black')\n",
    "    plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}'))\n",
    "    plt.savefig('figures/liver-plasma-corr-{}-{}.pdf'.format(str(genename) + '_' + protein, color_scheme), dpi=210, bbox_inches = 'tight') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data=data_corr.copy().set_index('Protein ID')\n",
    "to_plot = ['A0A286YEY1','P10643', 'P02741']\n",
    "for i in to_plot:\n",
    "    batch_plot_corr(i, df=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distribution of disease severity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In study participants whose plasma samples were analyzed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_scores = ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']\n",
    "biopsy_ids_mask = labtest['kleiner'].isin([0, 1, 2, 3, 4])\n",
    "biopsy_ids = labtest[biopsy_ids_mask]['Participant ID']\n",
    "plasma_cohortDist = {}\n",
    "for i in hist_scores:\n",
    "    plasma_cohortDist[i] = pd.DataFrame(labtest.set_index('Participant ID').reindex(biopsy_ids)[i].value_counts(1).round(2))\n",
    "    plasma_cohortDist[i] = plasma_cohortDist[i].sort_index()*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "for score in hist_scores[:1]:\n",
    "    df = plasma_cohortDist[score]\n",
    "    fig, ax = plt.subplots(figsize = (4, 4))\n",
    "    height = df[score]\n",
    "    plt.bar(x = np.arange(df.shape[0]), height = height, \n",
    "            facecolor = 'white', edgecolor = 'black', width = 0.5)\n",
    "    for i, v in enumerate(height):\n",
    "        ax.text(i-0.3 , v + 1, s = str(round(v, 1)) + '%', color='darkred', fontsize = 15)\n",
    "    plt.xticks(ticks = np.arange(df.shape[0]), fontsize = 15)\n",
    "    plt.yticks(fontsize = 15)\n",
    "    plt.ylim(0, height.max()*1.2)\n",
    "    plt.ylabel(score, fontsize = 15)\n",
    "    plt.title(score, fontsize = 15)\n",
    "    plt.savefig('figures/{}.pdf'.format(score), dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In study participants whose liver biopsy samples were analyzed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "liverTissue_PIDs = annotation_file[(annotation_file['Sample type'] == 'Liver') & (annotation_file['Group2'] == 'ALD')]['Participant ID'].to_list()\n",
    "\n",
    "hist_scores = ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']\n",
    "liver_cohortDist = {}\n",
    "for i in hist_scores:\n",
    "    liver_cohortDist[i] = pd.DataFrame(labtest.set_index('Participant ID').reindex(liverTissue_PIDs)[i].value_counts(1).round(2))\n",
    "    liver_cohortDist[i] = liver_cohortDist[i].sort_index()*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for score in hist_scores[:1]:\n",
    "    df = liver_cohortDist[score]\n",
    "    fig, ax = plt.subplots(figsize = (4, 4))\n",
    "    height = df[score]\n",
    "    plt.bar(x = np.arange(df.shape[0]), height = height, \n",
    "            facecolor = 'white', edgecolor = 'black', width = 0.5)\n",
    "    for i, v in enumerate(height):\n",
    "        ax.text(i-0.3 , v + 1, s = str(round(v, 0)) + '%', color='darkred', fontsize = 15)\n",
    "    plt.xticks(ticks = np.arange(df.shape[0]), fontsize = 15)\n",
    "    plt.yticks(fontsize = 15)\n",
    "    plt.ylim(0, height.max()*1.2)\n",
    "    plt.ylabel(score, fontsize = 15)\n",
    "    plt.title(score, fontsize = 15)\n",
    "    #plt.savefig('Clinical_parameters/tissueCohort/{}.pdf'.format(score), dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IDmapping_liverSampleID_to_ParticipantID = dict(zip(df_cli_liver.index, df_cli_liver['Participant ID']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Partial correlation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "## Prepare datasets\n",
    "- no imputation\n",
    "- no GALA-HP cohort or the at-risk group in ALD cohort\n",
    "- only study subjects who had histology examination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plasma_proteome = data_plasma_filtered_log_imputed.copy().T\n",
    "df_clinical = labtest.copy()\n",
    "df_clinical = df_clinical[np.logical_not(df_clinical.project=='GALA-HP')]\n",
    "col_to_keep = ['Participant ID', 'project', 'Sample ID', 'kleiner', 'nas_inflam', 'nas_steatosis_ordinal', 'age', 'bmi', 'gender_num', 'abstinent_num']\n",
    "df_clinical = df_clinical[col_to_keep].set_index('Sample ID')\n",
    "overlap_parID = intersection(df_plasma_proteome.index.to_list(), df_clinical.index.to_list())\n",
    "data_plasma_ptcorr = df_plasma_proteome.loc[overlap_parID].join(other = df_clinical)\n",
    "data_plasma_ptcorr = data_plasma_ptcorr[np.logical_not(data_plasma_ptcorr.kleiner==0.5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Liver dataset for partial correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_liver_proteome = data_liver_filtered_log_imputed.copy().T\n",
    "overlap_parID = intersection(df_liver_proteome.index.to_list(), df_cli_liver.index.to_list())\n",
    "col_to_keep = ['Participant ID', 'kleiner', 'nas_inflam', 'nas_steatosis_ordinal', 'age', 'bmi', 'gender_num', 'abstinent_num']\n",
    "data_liver_ptcorr = df_liver_proteome.loc[overlap_parID].join(other=df_cli_liver[col_to_keep])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": false
   },
   "source": [
    "## Compute partial correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "RECALCULATE_PARTIAL_CORR = False\n",
    "col_names = ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']\n",
    "FILE_RESULTS = os.path.join('data/processed/Partial_correlation_results.xlsx')\n",
    "\n",
    "if not RECALCULATE_PARTIAL_CORR:\n",
    "    ptcorr_liver = {}\n",
    "    ptcorr_plasma ={}\n",
    "    for i in col_names:\n",
    "        ptcorr_liver[i]=pd.read_excel(FILE_RESULTS, sheet_name='liver_'+i, index_col='protein')\n",
    "        ptcorr_plasma[i]=pd.read_excel(FILE_RESULTS, sheet_name='plasma_'+i, index_col='protein')\n",
    "        \n",
    "else:\n",
    "    data = data_liver_ptcorr.copy()\n",
    "    x_cols = data_liver_filtered.index.to_list()\n",
    "    np.random.seed(123)\n",
    "    ptcorr_liver = {}\n",
    "    for i in col_names:\n",
    "        if i!='nas_steatosis_ordinal':\n",
    "            ptcorr_liver[i] = partial_corr_pg(data=data, x_cols=x_cols, y_col=i, covar=covar)\n",
    "        elif i=='nas_steatosis_ordinal':\n",
    "            ptcorr_liver[i] = partial_corr_pg(data=data, x_cols=x_cols, y_col=i, covar=covar_steatosis)\n",
    "        \n",
    "    np.random.seed(123)\n",
    "    data = data_plasma_ptcorr.copy()\n",
    "    x_cols = data_plasma_filtered.index.to_list()\n",
    "    ptcorr_plasma = {}\n",
    "    for i in col_names:\n",
    "        if i!='nas_steatosis_ordinal':\n",
    "            ptcorr_plasma[i] = partial_corr_pg(data=data, x_cols=x_cols, y_col=i, covar=covar)\n",
    "        elif i=='nas_steatosis_ordinal':\n",
    "            ptcorr_plasma[i] = partial_corr_pg(data=data, x_cols=x_cols, y_col=i, covar=covar_steatosis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_RESULTS = os.path.join('data/processed/Partial_correlation_results.xlsx')\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    for i in col_names:\n",
    "        ptcorr_liver[i].to_excel(writer, sheet_name='liver_'+i)\n",
    "        ptcorr_plasma[i].to_excel(writer, sheet_name ='plasma_' + i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract significant correlations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Liver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_names = ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptcorr_liver_sig={}\n",
    "for i in col_names:\n",
    "    new_name = i+'_sig'\n",
    "    ptcorr_liver_sig[new_name] = return_partialcorr_sig(ptcorr_liver[i], thresh=0.3)\n",
    "    ptcorr_liver_sig[new_name]['Gene name'] = ptcorr_liver_sig[new_name].index.map(IDmapping_Liver_UniprotID_to_Genename)\n",
    "    ptcorr_liver_sig[new_name]['predicted to be secreted']=ptcorr_liver_sig[new_name]['Gene name'].map(secreted_dict)\n",
    "    ptcorr_liver_sig[new_name]['Protein IDs']=ptcorr_liver_sig[new_name].index.map(IDmapping_liver_ProteinID_to_ProteinIDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in col_names:\n",
    "    print(i + str(ptcorr_liver_sig[i+'_sig'].shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plasma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptcorr_plasma_sig={}\n",
    "for i in col_names:\n",
    "    new_name = i+'_sig'\n",
    "    ptcorr_plasma_sig[new_name] = return_partialcorr_sig(ptcorr_plasma[i], thresh=0.3)\n",
    "    ptcorr_plasma_sig[new_name]['Gene name'] = ptcorr_plasma_sig[new_name].index.map(IDmapping_UniprotID_to_Genename)\n",
    "    ptcorr_plasma_sig[new_name]['predicted to be secreted']=ptcorr_plasma_sig[new_name]['Gene name'].map(secreted_dict)\n",
    "    ptcorr_plasma_sig[new_name]['Protein IDs']=ptcorr_plasma_sig[new_name].index.map(IDmapping_plasma_ProteinID_to_ProteinIDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in col_names:\n",
    "    print(i + str(ptcorr_plasma_sig[i+'_sig'].shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_RESULTS = os.path.join('data/processed/Partial_correlation_sig.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    for i in ['kleiner_sig', 'nas_inflam_sig', 'nas_steatosis_ordinal_sig']:\n",
    "        ptcorr_plasma_sig[i].to_excel(writer, sheet_name='plasma' + i)\n",
    "        ptcorr_liver_sig[i].to_excel(writer, sheet_name='liver' + i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Top 20 plasma proteins correlated with steatosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=ptcorr_plasma['nas_steatosis_ordinal']\n",
    "df['Gene name']=df.index.map(IDmapping_UniprotID_to_Genename)\n",
    "df['predicted to be secreted']=df['Gene name'].map(secreted_dict)\n",
    "df['Protein IDs']=df.index.map(IDmapping_plasma_ProteinID_to_ProteinIDs)\n",
    "df=df.sort_values(by='r', ascending=False)\n",
    "df.to_csv('data/processed/partial_correlation_steatosis.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "- top 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sns.set_style('ticks')\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (9, 6))\n",
    "dfs = [ptcorr_plasma_sig[i][:20] for i in ['kleiner_sig', 'nas_inflam_sig', 'nas_steatosis_ordinal_sig']]\n",
    "axs = [ax1, ax2, ax3]\n",
    "titles = ['Fibrosis\\n(Kleiner)', 'Inflammation\\n(NAS_inflammation)', 'Steatosis\\n(NAS_steatosis)']\n",
    "for i in range(3):\n",
    "    ax = axs[i]\n",
    "    df = dfs[i]\n",
    "    color = np.where(df['predicted to be secreted']=='secreted', 'steelblue', 'white')\n",
    "    title = titles[i]\n",
    "    ax.barh(y = np.arange(20), width = df['r'][::-1], edgecolor = 'b', height = 0.65, color = color[::-1], hatch='/////')\n",
    "    ax.set_title(title, fontsize = 20)\n",
    "    ax.set_yticklabels(labels = df['Gene name'][::-1], fontsize = 13)\n",
    "    ax.set_xlabel('Spearman Corr. Coeff.', fontsize = 15)\n",
    "    ax.set_yticks(np.arange(20))\n",
    "    ax.set_xlim(0, 1)\n",
    "plt.subplots_adjust(wspace = 0.8)\n",
    "#fig.suptitle('Correlation to histology scores', fontsize = 25);\n",
    "plt.savefig('figures/PartialCorrelation/partial_correlation_plasma.pdf', dpi = 120, bbox_inches = 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- top 5 in box plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### export source data for figure 2e and figure S5e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top4_l_f = ptcorr_liver_sig['kleiner_sig'].sort_values(by='r', ascending=False).index[:4].tolist()\n",
    "top4_l_i = ptcorr_liver_sig['nas_inflam_sig'].sort_values(by='r', ascending=False).index[:4].tolist()\n",
    "top4_l_s  =ptcorr_liver_sig['nas_steatosis_ordinal_sig'].sort_values(by='r', ascending=False).index[:4].tolist()\n",
    "top_l = list(set(top4_l_f + top4_l_i + top4_l_s)) + ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']\n",
    "\n",
    "top4_p_f = ptcorr_plasma_sig['kleiner_sig'].sort_values(by='r', ascending=False).index[:4].tolist()\n",
    "top4_p_i = ptcorr_plasma_sig['nas_inflam_sig'].sort_values(by='r', ascending=False).index[:4].tolist()\n",
    "top4_p_s  =ptcorr_plasma_sig['nas_steatosis_ordinal_sig'].sort_values(by='r', ascending=False).index[:4].tolist()\n",
    "top_p = list(set(top4_p_f + top4_p_i + top4_p_s)) + ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sourcedata_fig2e = data_liver_ptcorr[top_l].melt(id_vars=['kleiner', 'nas_inflam', 'nas_steatosis_ordinal'], value_name='Log2 intensity', var_name='Protein ID')\n",
    "#de-couple histology scores\n",
    "sourcedata_fig2e = sourcedata_fig2e.melt(id_vars=['Protein ID', 'Log2 intensity'], var_name='histology', value_name='score')\n",
    "sourcedata_fig2e['Gene name']=sourcedata_fig2e['Protein ID'].map(IDmapping_Liver_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sourcedata_figs5e = data_plasma_ptcorr[top_p].melt(id_vars=['kleiner', 'nas_inflam', 'nas_steatosis_ordinal'], value_name='Log2 intensity', var_name='Protein ID')\n",
    "#de-couple histology scores\n",
    "sourcedata_figs5e = sourcedata_figs5e.melt(id_vars=['Protein ID', 'Log2 intensity'], var_name='histology', value_name='score')\n",
    "sourcedata_figs5e['Gene name']=sourcedata_figs5e['Protein ID'].map(IDmapping_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SOURCE_DATA_FIG2 = os.path.join('data/processed/source_data_figure2e_s5e.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(SOURCE_DATA_FIG2) as writer:\n",
    "    sourcedata_fig2e.to_excel(writer, sheet_name='Fig2_e')\n",
    "    sourcedata_figs5e.to_excel(writer, sheet_name='FigS5_e')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs = [ptcorr_plasma_sig[i] for i in ['kleiner_sig', 'nas_inflam_sig', 'nas_steatosis_ordinal_sig']]\n",
    "df1 = data_plasma_ptcorr.copy()\n",
    "df2 = dfs[2]\n",
    "\n",
    "top_df2 = df2[:4].index.to_list()\n",
    "parameter = 'nas_steatosis_ordinal'\n",
    "for protein in top_df2:\n",
    "    r = df2.loc[protein].r\n",
    "    #GN = IDmapping_Liver_UniprotID_to_Genename[protein]\n",
    "    GN = IDmapping_UniprotID_to_Genename[protein]\n",
    "    x = df1[parameter]\n",
    "    y = df1[protein]\n",
    "    fig, ax = plt.subplots(figsize = (3,3))\n",
    "    params = {'xtick.labelsize':15, \n",
    "             'ytick.labelsize':15,\n",
    "             'axes.labelsize':15,\n",
    "             'axes.titlesize':15}\n",
    "    sns.boxplot(x = x, y = y, color = 'white', width = 0.6, )\n",
    "    sns.stripplot(x=x, y=y, color='darkblue', size=3)\n",
    "    plt.title(GN, fontsize = 24)\n",
    "    plt.ylabel('MS signal [Log2]', fontsize = 20)\n",
    "    plt.yticks(fontsize=20)\n",
    "    plt.xticks(fontsize=20)\n",
    "    plt.annotate(s = 'R={}'.format(str(r.round(2))), xy = (0.05, 0.9), xycoords = 'axes fraction', fontsize = 18)\n",
    "    plt.rcParams.update(params)\n",
    "    plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}'))\n",
    "    plt.savefig('figures/PartialCorrelation/plasma_{}_{}_{}.pdf'.format(parameter, GN, protein), dpi = 120, bbox_inches= 'tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc-hr-collapsed": true
   },
   "source": [
    "# Global correlation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plasma = data_plasma_filtered_log_imputed.copy().rename_axis('Sample ID', axis = 1)\n",
    "df_plasma = df_plasma.rename(mapper = IDmapping_SampleID_to_ParticipantID, axis = 1)\n",
    "df_liver = data_liver_filtered_log_imputed.copy().rename_axis('Sample ID', axis = 1)\n",
    "df_liver = df_liver.rename(mapper = IDmapping_liverSampleID_to_ParticipantID, axis = 1)\n",
    "df_cli = labtest_numeric.copy()\n",
    "df_cli['Participant ID'] = df_cli.index.map(IDmapping_SampleID_to_ParticipantID)\n",
    "df_cli.set_index('Participant ID', inplace=True)\n",
    "\n",
    "data_corr_liver=df_liver.T.join(other=df_cli, how='inner')\n",
    "data_corr_plasma=df_plasma.T.join(other=df_cli, how='inner')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation of coagulation factors to INR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coag_factors = ['P00734', 'E7END6', 'P00742', 'P03951', 'K7ER74', 'B0YIW2', 'K7ERI9', 'inr']\n",
    "data_plasma_corr_inr = data_corr_plasma[coag_factors].dropna()\n",
    "plasma_corr_inr=pg.pairwise_corr(data=data_plasma_corr_inr, columns=['inr'])\n",
    "bins = [0.65, 0.85, 0.95, 1.05, 1.25, 2.7]\n",
    "data_plasma_corr_inr['inr_bins'] = pd.cut(data_plasma_corr_inr['inr'], bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Plotting bar plots of coagulation factors as a function of binned INR values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for proteinID in [i for i in coag_factors[:1] if i!='inr']: \n",
    "    GN = IDmapping_UniprotID_to_Genename[proteinID]\n",
    "    r_df = plasma_corr_inr.set_index('Y')\n",
    "    r = r_df.loc[proteinID]['r']\n",
    "    fig, ax = plt.subplots(figsize = (4,4))\n",
    "    sns.boxplot(x = data_plasma_corr_inr['inr_bins'], y = data_plasma_corr_inr[proteinID], color = 'royalblue', width = 0.6)\n",
    "    #plt.xlim(0.65, 1.65)\n",
    "    plt.ylabel(GN)\n",
    "    plt.xlabel('INR')\n",
    "    plt.xticks(rotation = 30)\n",
    "    plt.title('{} vs. INR'.format(GN))\n",
    "    plt.annotate(s = 'r = ' + str(r.round(2)), xy=[0.72, 0.92], xycoords = 'axes fraction', fontsize =16)\n",
    "    #plt.savefig('Clinical_parameters/figures/inr/{}_INR.png'.format(proteinID), dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "toc-autonumbering": true,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}