code/benchmarking/01.Compute-Binary-based-scores.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ff02550c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cell2cell as c2c\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from collections import defaultdict\n",
    "import random\n",
    "import itertools\n",
    "from tqdm.auto import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "888c8497",
   "metadata": {},
   "source": [
    "# Setups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "854c3c9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "rnaseq_setup = dict()\n",
    "rnaseq_setup['gene_col'] = 'gene_id'\n",
    "rnaseq_setup['drop_nangenes'] = True\n",
    "rnaseq_setup['log_transform'] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "368c984b",
   "metadata": {},
   "outputs": [],
   "source": [
    "meta_setup = dict()\n",
    "meta_setup['sample_col'] = '#SampleID'\n",
    "meta_setup['group_col'] = 'Groups'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "77dc0a5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ppi_setup = dict()\n",
    "ppi_setup['protein_cols'] = ['Ligand_WB', 'Receptor_WB']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1cf13724",
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoff_setup = dict()\n",
    "cutoff_setup['type'] = 'constant_value'\n",
    "cutoff_setup['parameter'] = 10 # TPM (Changed later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a6e9ea6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_setup = dict()\n",
    "analysis_setup['communication_score'] = 'expression_thresholding'\n",
    "analysis_setup['cci_score'] = 'count' # Changed later\n",
    "analysis_setup['cci_type'] = 'undirected'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac48582",
   "metadata": {},
   "source": [
    "# Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fe6b06e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_folder = '../../data/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b195f163",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "output_folder = '../../results/Benchmarking/'\n",
    "if not os.path.isdir(output_folder):\n",
    "    os.makedirs(output_folder, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "66026371",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = dict()\n",
    "files['rnaseq'] = data_folder + '/RNA-Seq/Celegans_RNASeqData_Cell.xlsx'\n",
    "files['metadata'] = data_folder + '/RNA-Seq/Celegans_cell_metadata.tsv'\n",
    "files['ppi'] = data_folder + '/PPI-Networks/Celegans-Curated-LR-pairs.xlsx'\n",
    "files['output_folder'] = output_folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "bad1515d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "if not os.path.isdir(files['output_folder']):\n",
    "    os.mkdir(files['output_folder'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e541a8fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening RNAseq datasets from /Users/earmingol/Dropbox/Universidad/UCSanDiego/Lab_Lewis/CCI-Project/Celegans/Project-cell2cell/Celegans-cell2cell-2022/Data//RNA-Seq/Celegans_RNASeqData_Cell.xlsx\n",
      "/Users/earmingol/Dropbox/Universidad/UCSanDiego/Lab_Lewis/CCI-Project/Celegans/Project-cell2cell/Celegans-cell2cell-2022/Data//RNA-Seq/Celegans_RNASeqData_Cell.xlsx was correctly loaded\n"
     ]
    }
   ],
   "source": [
    "rnaseq_data = c2c.io.load_rnaseq(rnaseq_file=files['rnaseq'],\n",
    "                                 gene_column=rnaseq_setup['gene_col'],\n",
    "                                 drop_nangenes=rnaseq_setup['drop_nangenes'],\n",
    "                                 log_transformation=rnaseq_setup['log_transform'],\n",
    "                                 format='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0c62ec59",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/earmingol/Dropbox/Universidad/UCSanDiego/Lab_Lewis/CCI-Project/Celegans/Project-cell2cell/Celegans-cell2cell-2022/Data//RNA-Seq/Celegans_cell_metadata.tsv was correctly loaded\n"
     ]
    }
   ],
   "source": [
    "meta = c2c.io.load_metadata(metadata_file=files['metadata'],\n",
    "                            cell_labels=rnaseq_data.columns,\n",
    "                            format='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "808ed2e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Opening PPI datasets from /Users/earmingol/Dropbox/Universidad/UCSanDiego/Lab_Lewis/CCI-Project/Celegans/Project-cell2cell/Celegans-cell2cell-2022/Data//PPI-Networks/Celegans-Curated-LR-pairs.xlsx\n",
      "/Users/earmingol/Dropbox/Universidad/UCSanDiego/Lab_Lewis/CCI-Project/Celegans/Project-cell2cell/Celegans-cell2cell-2022/Data//PPI-Networks/Celegans-Curated-LR-pairs.xlsx was correctly loaded\n",
      "Removing bidirectionality of PPI network\n",
      "Simplying PPI network\n"
     ]
    }
   ],
   "source": [
    "ppi_data = c2c.io.load_ppi(ppi_file=files['ppi'],\n",
    "                           interaction_columns=ppi_setup['protein_cols'],\n",
    "                           rnaseq_genes=list(rnaseq_data.index),\n",
    "                           format='auto')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "3629a661",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/earmingol/Dropbox/Universidad/UCSanDiego/Lab_Lewis/CCI-Project/Celegans/Project-cell2cell/Celegans-cell2cell-2022/Data//PPI-Networks/Celegans-Curated-LR-pairs.xlsx was correctly loaded\n"
     ]
    }
   ],
   "source": [
    "lr_pairs = c2c.io.load_table(files['ppi'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "b9de5eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_genes = list(set(lr_pairs['Ligand_WB'].unique()).union(set(lr_pairs['Receptor_WB'].unique())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "179b81d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "excluded_cells = ['Distal_tip_cells',\n",
    "                  'Sex_myoblasts',\n",
    "                  'Socket_cells',\n",
    "                  'Vulval_precursors',\n",
    "                  'flp-1(+)_interneurons']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "e1855a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "included_cells = list((set(rnaseq_data.columns)) - set(excluded_cells))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "52740d8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "rnaseq_data = rnaseq_data[included_cells]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bba7b5c9",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "ef09f303",
   "metadata": {},
   "outputs": [],
   "source": [
    "if analysis_setup['cci_type'] == 'undirected':\n",
    "    bi_ppi_data = c2c.preprocessing.bidirectional_ppi_for_cci(ppi_data=ppi_data, verbose=False)\n",
    "    ref_ppi = ppi_data\n",
    "else:\n",
    "    bi_ppi_data = ppi_data.copy()\n",
    "    ref_ppi = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "c9771407",
   "metadata": {},
   "outputs": [],
   "source": [
    "tpm_values = [5, 10, 50, 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "1458da26",
   "metadata": {},
   "outputs": [],
   "source": [
    "cci_scores = ['bray_curtis', 'count']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "52d47af3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Smillie setups\n",
    "permutations = 9999\n",
    "seed = 888"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "c2b9ec9d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cfa092322ca64609aad9bded03e9c44e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/9999 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "80d4b77fd5854ba2ab24152a7c20946e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/9999 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d5018a83a4564c9191a280cf0e39b117",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/9999 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be40e4e54fe54f578d6a8552101562d0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/9999 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for score in cci_scores:\n",
    "    analysis_setup['cci_score'] = score\n",
    "    for tpm in tpm_values:\n",
    "        cutoff_setup['parameter'] = tpm\n",
    "        \n",
    "        # Base scores\n",
    "        interaction_space = c2c.core.InteractionSpace(rnaseq_data=rnaseq_data[included_cells],\n",
    "                                                      ppi_data=bi_ppi_data,\n",
    "                                                      gene_cutoffs=cutoff_setup,\n",
    "                                                      communication_score=analysis_setup['communication_score'],\n",
    "                                                      cci_score=analysis_setup['cci_score'],\n",
    "                                                      cci_type=analysis_setup['cci_type'],\n",
    "                                                      verbose=False)\n",
    "\n",
    "        # Compute interactions\n",
    "        interaction_space.compute_pairwise_cci_scores(use_ppi_score=False, verbose=False)\n",
    "        cci_matrix = interaction_space.interaction_elements['cci_matrix'].copy()\n",
    "        \n",
    "        cci_matrix.to_csv(files['output_folder'] + 'CCI-Matrix-{}-{}-TPM.csv'.format(score.capitalize(), tpm))\n",
    "        \n",
    "        \n",
    "        ## Compute Smillie when score == 'count' ##\n",
    "        if score == 'count':\n",
    "            # Generate degree distribution for permutation\n",
    "            rnaseq_df = rnaseq_data.loc[lr_genes, :]\n",
    "            \n",
    "            cell_degrees = rnaseq_df.gt(cutoff_setup['parameter']).sum().sort_values()\n",
    "\n",
    "            min_val = rnaseq_df.gt(cutoff_setup['parameter']).sum().min()\n",
    "            max_val = rnaseq_df.gt(cutoff_setup['parameter']).sum().max()\n",
    "\n",
    "            degree_groups = defaultdict(set)\n",
    "            step = 10 # Bin size for same distribution\n",
    "\n",
    "            bins = np.arange(min_val, max_val, step)\n",
    "            assignment = np.digitize(cell_degrees.values, bins, right=False)\n",
    "\n",
    "            for i, k in enumerate(assignment):\n",
    "                degree_groups[k].add(list(cell_degrees.index)[i])\n",
    "\n",
    "            degree_groups = {k : list(v) for k, v in degree_groups.items()}\n",
    "\n",
    "            # Permutations\n",
    "            og_list = list(itertools.chain(*degree_groups.values()))\n",
    "            randomized_scores = []\n",
    "\n",
    "            for i in tqdm(range(permutations)):\n",
    "\n",
    "                shuffled_groups = []\n",
    "\n",
    "                for v in degree_groups.values():\n",
    "                    random.seed(seed+i)\n",
    "                    shuffled_groups.extend(random.sample(v, len(v)))\n",
    "\n",
    "                cell_mapper = dict(zip(og_list, shuffled_groups))\n",
    "\n",
    "                shuffled_rnaseq = rnaseq_df.copy()\n",
    "                shuffled_rnaseq = shuffled_rnaseq.rename(columns=cell_mapper)\n",
    "                shuffled_rnaseq = shuffled_rnaseq[included_cells] # Order of labels is important\n",
    "\n",
    "                # Compute CCI scores\n",
    "                interaction_space = c2c.core.InteractionSpace(rnaseq_data=shuffled_rnaseq,\n",
    "                                                          ppi_data=bi_ppi_data,\n",
    "                                                          gene_cutoffs=cutoff_setup,\n",
    "                                                          communication_score=analysis_setup['communication_score'],\n",
    "                                                          cci_score=analysis_setup['cci_score'],\n",
    "                                                          cci_type=analysis_setup['cci_type'],\n",
    "                                                          verbose=False)\n",
    "\n",
    "                # Compute interactions\n",
    "                interaction_space.compute_pairwise_cci_scores(use_ppi_score=False, verbose=False)\n",
    "\n",
    "                randomized_scores.append(interaction_space.interaction_elements['cci_matrix'].values.flatten())\n",
    "\n",
    "            # Add base scores as part of the null distribution\n",
    "            base_scores = cci_matrix.values.flatten()\n",
    "            randomized_scores.append(base_scores)\n",
    "            randomized_scores = np.array(randomized_scores)\n",
    "\n",
    "            # Compute P-values\n",
    "            pvals = np.ones(base_scores.shape)\n",
    "\n",
    "            for i in range(len(base_scores)):\n",
    "                dist = randomized_scores[:, i]\n",
    "                pvals[i] = c2c.stats.permutation.compute_pvalue_from_dist(obs_value=base_scores[i],\n",
    "                                                                dist=dist,\n",
    "                                                                consider_size=True,\n",
    "                                                                comparison='upper'\n",
    "                                                                )\n",
    "            pval_df = pd.DataFrame(pvals.reshape(cci_matrix.shape), index=cci_matrix.index, \n",
    "                                   columns=cci_matrix.columns)\n",
    "\n",
    "            output_df = -np.log10(pval_df)\n",
    "            \n",
    "            output_df.to_csv(files['output_folder'] + 'CCI-Matrix-Smillie-{}-TPM.csv'.format(tpm))\n",
    "            output_df.to_csv(data_folder + '/Benchmarking/CCI-Matrix-Smillie-{}-TPM.csv'.format(tpm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84d20e24",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}