analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Statistical analysis\n",
    "\n",
    "_Michał Denkiewicz 2024_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import string\n",
    "import functools\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import scipy.stats\n",
    "import statsmodels.stats.multitest\n",
    "import networkx as nx\n",
    "import pyranges as pr\n",
    "\n",
    "from itertools import product\n",
    "from functools import partial\n",
    "\n",
    "from knots_tools import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Number of minors detected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR = './data'\n",
    "\n",
    "ccds = pd.read_csv(\n",
    "     os.path.join(DATA_DIR, 'all_ccds.csv'),\n",
    "     index_col=CCD_INDEX_NAMES,\n",
    "     dtype=CCD_INDEX_DTYPES\n",
    ")\n",
    "ccds.info()\n",
    "ccds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minors = pd.read_csv(\n",
    "    os.path.join(DATA_DIR, 'all_minors.csv'),\n",
    "    index_col=CCD_INDEX_NAMES + ['minor_id'],\n",
    "    dtype=CCD_INDEX_DTYPES\n",
    ")\n",
    "minors.info()\n",
    "minors.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minors.groupby(['dataset'], observed=True).size()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minors.groupby(['dataset', 'chromosome'], observed=True).size().unstack('dataset', fill_value=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def pretty_formated_percent(x, precision=1):\n",
    "    m = x.mean()\n",
    "    if np.isnan(m):\n",
    "        m = 0.0\n",
    "    return f'{100.0 * m:.{precision}f}%'\n",
    "\n",
    "ccds.has_minors.groupby(['dataset'], observed=True).agg(['sum', pretty_formated_percent]).rename(columns={'sum': 'N', 'pretty_formated_percent': 'percent'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_tab = ccds.has_minors.groupby(['dataset', 'chromosome'], observed=True).agg(['sum', pretty_formated_percent]).rename(columns={'sum': 'N', 'pretty_formated_percent': 'percent'})\n",
    "_tab.columns.names = ['var']\n",
    "_tab = _tab.unstack(['dataset'], fill_value=0)\n",
    "_tab = _tab.reorder_levels(['dataset', 'var'], axis=1).sort_index(axis=1)\n",
    "_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Node / edge characteristics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes = pd.read_csv(\n",
    "    os.path.join(DATA_DIR, 'all_nodes.csv'),\n",
    "    index_col=CCD_INDEX_NAMES + ['node_id'],\n",
    "    dtype=CCD_INDEX_DTYPES\n",
    ")\n",
    "nodes.info()\n",
    "nodes.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "edges = pd.read_csv(\n",
    "    os.path.join(DATA_DIR, 'all_edges.csv'),\n",
    "    index_col=CCD_INDEX_NAMES + ['edge_id'],\n",
    "    dtype=CCD_INDEX_DTYPES\n",
    ")\n",
    "edges.info()\n",
    "edges.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes.in_minor.groupby(['dataset'], observed=True).agg(['sum', pretty_formated_percent, 'count']).rename(columns={'sum': 'N', 'pretty_formated_percent': 'percent', 'count': 'total'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "edges.in_minor.groupby(['dataset'], observed=True).agg(['sum', pretty_formated_percent, 'count']).rename(columns={'sum': 'N', 'pretty_formated_percent': 'percent', 'count': 'total'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graph measures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ccds.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "measure_labels = {\n",
    "    'n_nodes': 'No. Nodes',\n",
    "    'n_edges': 'No. Edges',\n",
    "    'density': 'Density',\n",
    "    'degree_mean': 'Avg. Degree',\n",
    "    'closeness_mean': 'Avg. Closeness',\n",
    "    'betweenness_mean': 'Avg. Betweenness'\n",
    "}\n",
    "\n",
    "data_long = ccds.set_index('has_minors', append=True).loc[:, list(measure_labels.keys())].stack().to_frame('value')\n",
    "data_long.index.names = data_long.index.names[:-1] + ['measure']\n",
    "data_long"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def do_test(df, var, by, which='med'):\n",
    "    (key0, g0), (key1, g1) = list(df[var].groupby(df[by]))\n",
    "    assert key0 == False and key1 == True\n",
    "    assert which.lower() in ('med', 't', 'u')\n",
    "    if which.lower() != 't':\n",
    "        m0 = g0.median()\n",
    "        m1 = g1.median()\n",
    "        mg = df[var].median()\n",
    "    else:\n",
    "        m0 = g0.mean()\n",
    "        m1 = g1.mean()\n",
    "        mg = df[var].mean()\n",
    "    try:\n",
    "        match which.lower():\n",
    "            case 'med':\n",
    "                stat, p, *_ = scipy.stats.median_test(g0, g1)\n",
    "            case 't':\n",
    "                stat, p, *_ = scipy.stats.ttest_ind(g0, g1, equal_var=False)\n",
    "            case 'u':\n",
    "                stat, p, *_ = scipy.stats.mannwhitneyu(g0, g1, use_continuity=True, alternative='two-sided')\n",
    "    except Exception as ex:\n",
    "        stat, p, mg = np.nan, np.nan, mg\n",
    "    return {'stat': stat, 'pvalue': p, 'm0': m0, 'm1': m1, 'gm': mg, 'n0': len(g0), 'n1': len(g1)}\n",
    "\n",
    "stats_tab = pd.DataFrame.from_dict(\n",
    "    data_long.reset_index().groupby(['dataset', 'measure'], observed=False).apply(do_test, 'value', 'has_minors', include_groups=False).to_dict(),\n",
    "    orient='index',    \n",
    ")\n",
    "stats_tab.index.names = ['dataset', 'measure']\n",
    "stats_tab['dir'] = np.where(stats_tab['m0'] > stats_tab['m1'], '>', '<')\n",
    "_, stats_tab['pcorr'], *_ = statsmodels.stats.multitest.multipletests(stats_tab['pvalue'], method='holm')\n",
    "stats_tab['sig'] = pd.cut(stats_tab['pcorr'], bins=[-1, 0.001, 0.01, 0.05, 1], labels=['***', '**', '*', 'ns'], right=True)\n",
    "stats_tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pretty_stats_tab = stats_tab.reset_index().replace({'measure': measure_labels, 'dataset': {'GM12878': 'GM12878(IS)', 'GM12878lr': 'GM12878(LR)'}}).set_index(['dataset', 'measure'])\n",
    "pretty_stats_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiplexes (ChIA-Drop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "fragments_df  = pd.read_csv(os.path.join(DATA_DIR, 'chiadrop_fragments.csv'))\n",
    "segments_df = pd.read_csv(os.path.join(DATA_DIR, 'all_minor_segments.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fragments_pr = pr.PyRanges(fragments_df.rename(columns={'chromosome': 'Chromosome', 'frag_start': 'Start', 'frag_end': 'End'}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "ccds_gm12878 = ccds.loc['GM12878']\n",
    "ccds_pr = pr.PyRanges(ccds_gm12878.reset_index().loc[:, ['ccd_id', 'chromosome', 'start', 'end', 'has_minors', 'length']].rename(columns={'chromosome': 'Chromosome', 'start': 'Start', 'end': 'End'}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "segments_gm12878 = segments_df.loc[segments_df.dataset == 'GM12878', :]\n",
    "segments_pr = pr.PyRanges(segments_gm12878.rename(columns={'chromosome': 'Chromosome', 'start': 'Start', 'end': 'End'}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ccd2frag = ccds_pr.join(fragments_pr, suffix='_frag', how='left').df\n",
    "ccd2frag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ccd2frag_stats = pd.concat({\n",
    "    'all': ccds_gm12878.groupby('has_minors').size(),\n",
    "    'has_fragments': ccd2frag.groupby(['has_minors']).ccd_id.nunique()\n",
    "}).unstack('has_minors', fill_value=0)\n",
    "ccd2frag_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy.stats.chi2_contingency(\n",
    "   ccd2frag_stats, correction=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### GEMs per unit length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gem_cnt = ccd2frag.groupby(['has_minors'], observed=True).gem_id.nunique()\n",
    "gem_cnt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tot_ccd_len = ccds_gm12878.set_index('has_minors', append=True).groupby('has_minors').length.sum()\n",
    "tot_ccd_len"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gems_per_mb = gem_cnt / (tot_ccd_len / 1e6)  # per 1Mb\n",
    "print(gems_per_mb[1] / gems_per_mb[0])\n",
    "gems_per_mb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Order of interactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arity_stats = ccd2frag.groupby(['n_fragments', 'has_minors']).gem_id.nunique().unstack('has_minors', fill_value=0)\n",
    "arity_stats.loc[10, :] += arity_stats.loc[11:, :].sum(axis=0)\n",
    "arity_stats = arity_stats.loc[:10, :]\n",
    "arity_stats['prop_has_minors'] = arity_stats[True] / (arity_stats[True] + arity_stats[False])\n",
    "arity_stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linked_ccds_multiplex_arity = ccd2frag.loc[ccd2frag.has_minors, :]\\\n",
    "    .groupby(['ccd_id', 'gem_id']).size().to_frame('n_mapped')\\\n",
    "    .groupby(['ccd_id', 'n_mapped']).size()\\\n",
    "    .unstack('n_mapped', fill_value=0)\n",
    "linked_ccds_multiplex_arity.loc[:, 6] += linked_ccds_multiplex_arity.loc[:, 7:].sum(axis=1)\n",
    "linked_ccds_multiplex_arity = linked_ccds_multiplex_arity.loc[:, :6]\n",
    "linked_ccds_multiplex_arity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(linked_ccds_multiplex_arity.loc[:, 6] > 0).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seg2frag = segments_pr.join(fragments_pr, suffix='_frag', how='left').df\n",
    "seg2frag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapped_segments = seg2frag.groupby(['Chromosome', 'ccd_id', 'minor_id', 'gem_id'], observed=True).segment_idx.nunique().to_frame('n_mapped').reset_index()\n",
    "mapped_segments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapped_segments.n_mapped.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_stat = mapped_segments.groupby(['Chromosome', 'ccd_id'], observed=True).n_mapped.max().value_counts().sort_index().to_frame('n_ccds')\n",
    "map_stat['p_ccds'] = np.round(100.0 * map_stat['n_ccds'] / len(mapped_segments.groupby(['Chromosome', 'ccd_id'], observed=True)), 1)\n",
    "print(np.cumsum(map_stat.p_ccds.to_numpy()[::-1]))\n",
    "map_stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fully_supported_minors = mapped_segments[mapped_segments['n_mapped'] >= 6]\n",
    "fully_supported_minors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ccds_gm12878.reset_index().rename(columns={'chromosome': 'Chromosome'})[['Chromosome', 'ccd_id', 'length', 'n_nodes', 'n_edges']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.merge(\n",
    "    fully_supported_minors,\n",
    "    ccds_gm12878.reset_index().rename(columns={'chromosome': 'Chromosome'})[['Chromosome', 'ccd_id', 'length', 'n_nodes', 'n_edges', 'start', 'end']],\n",
    "    on=['Chromosome', 'ccd_id'],\n",
    "    how='left'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Compartments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "raw_compartments = pd.read_csv('./raw_data/hg38_GM12878_compartments_250kb.bedgraph', header=None, sep='\\t')\n",
    "raw_compartments.columns = ['chromosome', 'start', 'end', 'score']\n",
    "raw_compartments['chromosome'] = raw_compartments['chromosome'].astype(HumanChromosomeDtype)\n",
    "raw_compartments = raw_compartments.loc[~raw_compartments.score.isna(), :]\n",
    "raw_compartments['compartment'] = np.where(raw_compartments.score >= 0, 'A', 'B')\n",
    "print(raw_compartments.groupby('compartment').size())\n",
    "\n",
    "comp_pr = pr.PyRanges(raw_compartments.rename(columns={'chromosome': 'Chromosome', 'start': 'Start', 'end': 'End'}))\n",
    "comp_pr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_compartments = pd.read_csv(\n",
    "    './raw_data/GSE63525_GM12878_subcompartments.bed',\n",
    "    header=None, sep='\\s+', usecols=list(range(4))\n",
    ").sort_values([0, 1, 2])\n",
    "raw_compartments.columns = ['chromosome', 'start', 'end', 'subcompartment']\n",
    "raw_compartments['end'] -= 1  # to prevent spurious overlaps\n",
    "raw_compartments['chromosome'] = raw_compartments['chromosome'].astype(HumanChromosomeDtype)\n",
    "raw_compartments['compartment'] = raw_compartments['subcompartment'].str[0]\n",
    "raw_compartments = raw_compartments.loc[~raw_compartments.compartment.isna(), :]\n",
    "\n",
    "print(raw_compartments.groupby('compartment').size())\n",
    "print(raw_compartments.groupby('subcompartment').size())\n",
    "\n",
    "comp_pr = pr.PyRanges(raw_compartments.rename(columns={'chromosome': 'Chromosome', 'start': 'Start', 'end': 'End'}))\n",
    "comp_pr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ccd2comp = ccds_pr.join(comp_pr, suffix='_comp', report_overlap=True).df.set_index(['Chromosome', 'ccd_id'])\n",
    "ccd2comp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_overlap_by_comp = ccd2comp.groupby(['has_minors', 'Chromosome', 'ccd_id', 'compartment'], observed=False)\\\n",
    "    .Overlap.sum().unstack('compartment', fill_value=0)\n",
    "coverage = ccds_gm12878.reset_index()\\\n",
    "    .rename(columns={'chromosome': 'Chromosome'})\\\n",
    "    .set_index(['has_minors', 'Chromosome', 'ccd_id'])\\\n",
    "    .loc[:, ['length']]\\\n",
    "    .join(_overlap_by_comp)\n",
    "coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coverage['A'] /= coverage['length']\n",
    "coverage['B'] /= coverage['length']\n",
    "coverage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coverage.groupby('has_minors').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coverage['compartment'] = 'mixed'\n",
    "coverage.loc[coverage['A'] >= 0.6, 'compartment'] = 'A'\n",
    "coverage.loc[coverage['B'] >= 0.6, 'compartment'] = 'B'\n",
    "tab_comp = coverage.groupby('has_minors').compartment.value_counts().to_frame('n_ccds')\n",
    "tab_comp['p_ccds'] = 100.0 * tab_comp['n_ccds'] / tab_comp.groupby('has_minors').n_ccds.sum()\n",
    "tab_comp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_cont_table = tab_comp.n_ccds.unstack('compartment', fill_value=0)\n",
    "print(_cont_table)\n",
    "scipy.stats.chi2_contingency(_cont_table, correction=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = tab_comp.p_ccds\n",
    "print(p.loc[False, 'B'] / p.loc[False, 'A'])\n",
    "print(p.loc[True, 'B'] / p.loc[True, 'A'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coverage['has_boundary'] = (coverage['A'] > 0.0) & (coverage['B'] > 0.0)\n",
    "tab_bound = coverage.groupby('has_minors').has_boundary.value_counts().to_frame('n_ccds')\n",
    "tab_bound['p_ccds'] = 100.0 * tab_bound['n_ccds'] / tab_bound.groupby('has_minors').n_ccds.sum()\n",
    "tab_bound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_cont_table = tab_bound.n_ccds.unstack('has_boundary', fill_value=0)\n",
    "print(_cont_table)\n",
    "scipy.stats.chi2_contingency(_cont_table, correction=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_dist_linked_to_nonlinked(X, var, xlab, yticks=False, palette='colorblind', scatter=False, ax=None, xlim=(0, None), log=False, legend=False):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    else:\n",
    "        plt.sca(ax)\n",
    "    hue_order = ['non-linked', 'linked']\n",
    "    sns.boxplot(\n",
    "        X, x=var, y=\"dataset\", hue=\"has_minors_str\",\n",
    "        whis=[0, 100], palette=palette, width=0.8,\n",
    "        ax=ax, linewidth=0.75, log_scale=log, hue_order=hue_order,\n",
    "        legend=legend\n",
    "    )\n",
    "    if scatter:\n",
    "        ptpal = sns.color_palette(['white', 'white'], as_cmap=True)\n",
    "        # ptpal = sns.color_palette(['gray', 'gray'], as_cmap=True)\n",
    "        # ptpal = 'tab10'\n",
    "        sns.stripplot(\n",
    "            X, x=var, y=\"dataset\", hue=\"has_minors_str\",\n",
    "            size=1, alpha=.2, palette=ptpal, dodge=True,\n",
    "            linewidth=0.25, edgecolor='auto', log_scale=log,\n",
    "            ax=ax, legend=False, hue_order=hue_order\n",
    "        )\n",
    "    if xlim is not None:\n",
    "        ax.set_xlim(xlim)\n",
    "    ax.xaxis.grid(True)\n",
    "    ax.set(xlabel=xlab)\n",
    "    ax.set(ylabel=\"\")\n",
    "    if not yticks:\n",
    "        plt.tick_params(labelleft=False, left=False)\n",
    "    if legend:\n",
    "        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0., title='')\n",
    "    sns.despine(trim=False, left=True)    \n",
    "\n",
    "ccd_stats = ccds.reset_index()\n",
    "ccd_stats['dataset_str'] = ccd_stats.dataset.astype('str')\n",
    "ccd_stats['has_minors_str'] = ccd_stats.has_minors.replace({True: 'linked', False: 'non-linked'})\n",
    "print(ccd_stats.dataset.value_counts())\n",
    "\n",
    "_plt = functools.partial(compare_dist_linked_to_nonlinked, ccd_stats, scatter=False)\n",
    "\n",
    "def _lab(s, x=-0.1, y=1.15):\n",
    "    ax = plt.gca()\n",
    "    ax.text(x, y, s, horizontalalignment='center', verticalalignment='center', transform=ax.transAxes, fontsize=14)\n",
    "\n",
    "matplotlib.rcParams.update({'font.size': 8})\n",
    "\n",
    "fig, axs = plt.subplots(\n",
    "    nrows=2, ncols=4, figsize=(7.8, 3.8), dpi=300,\n",
    "    sharex=False, sharey=False\n",
    ")\n",
    "iax = iter(axs.flatten())\n",
    "ilab = iter(string.ascii_uppercase)\n",
    "\n",
    "# A-F\n",
    "_plt('n_nodes', 'No. Nodes', log=True, xlim=(1, None), ax=next(iax), yticks=True); _lab(next(ilab), -.5)\n",
    "_plt('n_edges', 'No. Edges', log=True, xlim=(1, None), ax=next(iax)); _lab(next(ilab))\n",
    "_plt('density', 'Density', log=True, xlim=(1e-5, 1), ax=next(iax)); _lab(next(ilab))\n",
    "_plt('degree_mean', 'Avg. Degree', ax=next(iax), legend=True); _lab(next(ilab))\n",
    "_plt('closeness_mean', 'Avg. Closeness', ax=next(iax), yticks=True); _lab(next(ilab), -.5)\n",
    "_plt('betweenness_mean', 'Avg. Betweenness', ax=next(iax)); _lab(next(ilab))\n",
    "\n",
    "# G\n",
    "plt.sca(next(iax))\n",
    "sns.barplot(\n",
    "    x=map_stat.index, y='n_ccds', data=map_stat.reset_index(),\n",
    "    color='tan', linewidth=.5, edgecolor=\"black\"\n",
    ")\n",
    "_lab(next(ilab))\n",
    "plt.xlabel('Highest multiplex\\ninteraction order')\n",
    "plt.ylabel('No. of CCDs')\n",
    "\n",
    "# H\n",
    "plt.sca(next(iax))\n",
    "tmp = tab_comp.reset_index()\n",
    "tmp['has_minors'] = tmp['has_minors'].replace({False: 'non-linked', True: 'linked'})\n",
    "tmp = tmp.set_index(['has_minors', 'compartment'])\n",
    "tmp.p_ccds.unstack('compartment').plot(\n",
    "    kind='bar', stacked=True, ax=plt.gca(), color=sns.color_palette('Set2', 3),\n",
    "    linewidth=.5, edgecolor='black'\n",
    ")\n",
    "_lab(next(ilab))\n",
    "plt.ylim(0, 100)\n",
    "plt.xlabel('')\n",
    "plt.xticks(rotation=0)\n",
    "plt.ylabel('% CCDs')\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "# next(iax).axis('off')\n",
    "plt.subplots_adjust(hspace=0.7, wspace=0.5)\n",
    "\n",
    "plt.savefig(f'./plots/Figure_4_plots.svg', bbox_inches='tight')\n",
    "# plt.savefig(f'./plots/Figure_4_plots.png', bbox_inches='tight')\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sci1",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}