ETL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import requests\n",
    "import re\n",
    "from io import StringIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/luke/projects/cell-type-deconvolution/.venv/lib/python3.11/site-packages/openpyxl/worksheet/_reader.py:329: UserWarning: Unknown extension is not supported and will be removed\n",
      "  warn(msg)\n"
     ]
    }
   ],
   "source": [
    "signature_matrix_url = 'https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4739640/bin/NIHMS670442-supplement-2.xlsx'\n",
    "filename = 'signature_matrix'\n",
    "try:\n",
    "    open(f'./data/raw/{filename}.xlsx', 'r')\n",
    "except FileNotFoundError:\n",
    "    r = requests.get(signature_matrix_url)\n",
    "    open(f'./data/raw/{filename}.xlsx', 'wb').write(r.content)\n",
    "df = pd.read_excel(f'./data/raw/{filename}.xlsx', sheet_name='SuppTable1_GEP_Matrix', header=13, index_col='Gene symbol')\n",
    "df.drop(['Entrez ID', 'Max probeseta'], axis=1).to_csv(f'./data/{filename}.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/luke/projects/cell-type-deconvolution/.venv/lib/python3.11/site-packages/openpyxl/worksheet/_reader.py:329: UserWarning: Unknown extension is not supported and will be removed\n",
      "  warn(msg)\n"
     ]
    }
   ],
   "source": [
    "FILENAME = 'ground_truth_proportions' \n",
    "try:\n",
    "    open(f'./data/raw/{FILENAME}.xlsx', 'r')\n",
    "except FileNotFoundError:\n",
    "    r = requests.get('https://static-content.springer.com/esm/art%3A10.1038%2Fnmeth.3337/MediaObjects/41592_2015_BFnmeth3337_MOESM205_ESM.xlsx')\n",
    "    open(f'./data/raw/{FILENAME}.xlsx', 'wb').write(r.content)\n",
    "pd.read_excel(f'./data/raw/{FILENAME}.xlsx').to_csv(f'./data/{FILENAME}.csv', index=False)\n",
    "\n",
    "signature_matrix = pd.read_csv(f\"./data/signature_matrix.csv\", index_col=0)\n",
    "ground_truth = pd.read_csv('./data/ground_truth_proportions.csv').pivot(index='Sample ID', columns='Cell type', values='Flow (%)')\n",
    "signature_matrix.columns\n",
    "COL_MAPPING_GT_TO_SIG = {\n",
    "    'Activated memory CD4 T cells': 'T cells CD4 memory activated', \n",
    "    'CD8 T cells': 'T cells CD8',\n",
    "    'Memory B cells': 'B cells memory',\n",
    "    'Monocytes': 'Monocytes',\n",
    "    'NK cells': 'NK cells activated', \n",
    "    'Naïve B cells': 'B cells naive', \n",
    "    'Naïve CD4 T cells': 'T cells CD4 naive',\n",
    "    'Resting memory CD4 T cells': 'T cells CD4 memory resting',\n",
    "\n",
    "}\n",
    "new_ground_truth_cols = [COL_MAPPING_GT_TO_SIG[c] for c in ground_truth.columns]\n",
    "ground_truth.columns = new_ground_truth_cols\n",
    "assert sum(ground_truth.columns.isin(signature_matrix.columns) == False ) == 0\n",
    "ground_truth.to_csv('./data/ground_truth_proportions.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    file = open('./data/raw/HumanHT-12 V4.0 expression beadchip.txt', 'r').read()\n",
    "except FileNotFoundError:\n",
    "    r = requests.get('https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?mode=raw&is_datatable=true&acc=GPL10558&id=50081&db=GeoDb_blob135')\n",
    "    open('./data/raw/HumanHT-12 V4.0 expression beadchip.txt', 'wb').write(r.content)\n",
    "    file = r.text\n",
    "raw_data = file.split(\"#GB_ACC = \\n\")[1]\n",
    "beadchip_probe_definitions = pd.read_csv(StringIO(raw_data), sep='\\t', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I was also given this probe definition file, but I don't know its provenance\n",
    "alt_beadchip_probe_definitions = pd.read_csv('./data/HumanHT-12_v4.txt', index_col=0)\n",
    "alt_beadchip_probe_definitions = (alt_beadchip_probe_definitions\n",
    "    .reset_index()\n",
    "    .dropna()\n",
    "    .set_index('ILLUMINA HumanHT 12 V4 probe')\n",
    "    .groupby('ILLUMINA HumanHT 12 V4 probe').head(1)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_mix_ids = {\n",
    "    \"GSM1587800\": \"17-002\",\n",
    "    \"GSM1587801\": \"17-006\",\n",
    "    \"GSM1587802\": \"17-019\",\n",
    "    \"GSM1587803\": \"17-023\",\n",
    "    \"GSM1587804\": \"17-026\",\n",
    "    \"GSM1587805\": \"17-027\",\n",
    "    \"GSM1587806\": \"17-030\",\n",
    "    \"GSM1587807\": \"17-034\",\n",
    "    \"GSM1587808\": \"17-040\",\n",
    "    \"GSM1587809\": \"17-041\",\n",
    "    \"GSM1587810\": \"17-042\",\n",
    "    \"GSM1587811\": \"17-043\",\n",
    "    \"GSM1587812\": \"17-045\",\n",
    "    \"GSM1587813\": \"17-047\",\n",
    "    \"GSM1587814\": \"17-054\",\n",
    "    \"GSM1587815\": \"17-055\",\n",
    "    \"GSM1587816\": \"17-057\",\n",
    "    \"GSM1587817\": \"17-058\",\n",
    "    \"GSM1587818\": \"17-060\",\n",
    "    \"GSM1587819\": \"17-061\",\n",
    "}\n",
    "\n",
    "dfs = []\n",
    "for mix_id in cell_mix_ids.keys():\n",
    "    try:\n",
    "        html = open(f'./data/raw/cell_mixes/{mix_id}.html', 'r').read()\n",
    "    except FileNotFoundError:\n",
    "        html = requests.get(f'https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc={mix_id}').text\n",
    "        open(f'./data/raw/cell_mixes/{mix_id}.html', 'w').write(html)\n",
    "    split = html.split('#Detection Pval = \\n')[1].split('\\n')\n",
    "    assert len(split) == 47_336\n",
    "    raw_cols = split[0].split('</strong><strong>\\t')\n",
    "    clean = re.compile('<.*?>')\n",
    "    cols = '\\t'.join([re.sub(clean, '', string) for string in raw_cols])\n",
    "    df = pd.read_csv(StringIO('\\n'.join([cols] + split[2:47_326])), sep='\\t', dtype=str, index_col=0)\n",
    "    mix_name = cell_mix_ids[mix_id]\n",
    "    df = df.rename(columns={'VALUE': mix_name})\n",
    "    dfs.append(df[mix_name])\n",
    "all_mixes = pd.concat(dfs, axis=1)\n",
    "\n",
    "def get_HUGO_from_probe(probe):\n",
    "    try:\n",
    "        return alt_beadchip_probe_definitions.loc[probe]['HGNC symbol']\n",
    "        # return beadchip_probe_definitions.loc[probe]['Symbol']\n",
    "    except KeyError:\n",
    "        return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_mixes = all_mixes.astype(float)\n",
    "# Map to gene names from probe names using beadchip probe definitions\n",
    "all_mixes['Gene symbol'] = all_mixes.index.map(get_HUGO_from_probe)\n",
    "all_mixes.dropna(inplace=True, subset='Gene symbol')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>17-002</th>\n",
       "      <th>17-006</th>\n",
       "      <th>17-019</th>\n",
       "      <th>17-023</th>\n",
       "      <th>17-026</th>\n",
       "      <th>17-027</th>\n",
       "      <th>17-030</th>\n",
       "      <th>17-034</th>\n",
       "      <th>17-040</th>\n",
       "      <th>17-041</th>\n",
       "      <th>17-042</th>\n",
       "      <th>17-043</th>\n",
       "      <th>17-045</th>\n",
       "      <th>17-047</th>\n",
       "      <th>17-054</th>\n",
       "      <th>17-055</th>\n",
       "      <th>17-057</th>\n",
       "      <th>17-058</th>\n",
       "      <th>17-060</th>\n",
       "      <th>17-061</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ID_REF</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ILMN_1653618</th>\n",
       "      <td>85.132084</td>\n",
       "      <td>62.618970</td>\n",
       "      <td>63.335824</td>\n",
       "      <td>81.755322</td>\n",
       "      <td>68.364577</td>\n",
       "      <td>68.421480</td>\n",
       "      <td>76.784118</td>\n",
       "      <td>64.457203</td>\n",
       "      <td>102.787560</td>\n",
       "      <td>100.798411</td>\n",
       "      <td>83.929387</td>\n",
       "      <td>129.423574</td>\n",
       "      <td>108.943301</td>\n",
       "      <td>75.223829</td>\n",
       "      <td>83.607891</td>\n",
       "      <td>83.055035</td>\n",
       "      <td>99.582777</td>\n",
       "      <td>85.947656</td>\n",
       "      <td>90.397115</td>\n",
       "      <td>60.040611</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ILMN_2137536</th>\n",
       "      <td>42.611526</td>\n",
       "      <td>49.924148</td>\n",
       "      <td>48.327468</td>\n",
       "      <td>61.781762</td>\n",
       "      <td>40.804099</td>\n",
       "      <td>57.146188</td>\n",
       "      <td>50.908165</td>\n",
       "      <td>65.051624</td>\n",
       "      <td>60.183649</td>\n",
       "      <td>87.408230</td>\n",
       "      <td>46.616800</td>\n",
       "      <td>78.086649</td>\n",
       "      <td>124.652993</td>\n",
       "      <td>76.715498</td>\n",
       "      <td>67.866845</td>\n",
       "      <td>52.108945</td>\n",
       "      <td>58.353998</td>\n",
       "      <td>84.210700</td>\n",
       "      <td>63.367195</td>\n",
       "      <td>51.248614</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 17-002     17-006     17-019     17-023     17-026  \\\n",
       "ID_REF                                                                \n",
       "ILMN_1653618  85.132084  62.618970  63.335824  81.755322  68.364577   \n",
       "ILMN_2137536  42.611526  49.924148  48.327468  61.781762  40.804099   \n",
       "\n",
       "                 17-027     17-030     17-034      17-040      17-041  \\\n",
       "ID_REF                                                                  \n",
       "ILMN_1653618  68.421480  76.784118  64.457203  102.787560  100.798411   \n",
       "ILMN_2137536  57.146188  50.908165  65.051624   60.183649   87.408230   \n",
       "\n",
       "                 17-042      17-043      17-045     17-047     17-054  \\\n",
       "ID_REF                                                                  \n",
       "ILMN_1653618  83.929387  129.423574  108.943301  75.223829  83.607891   \n",
       "ILMN_2137536  46.616800   78.086649  124.652993  76.715498  67.866845   \n",
       "\n",
       "                 17-055     17-057     17-058     17-060     17-061  \n",
       "ID_REF                                                               \n",
       "ILMN_1653618  83.055035  99.582777  85.947656  90.397115  60.040611  \n",
       "ILMN_2137536  52.108945  58.353998  84.210700  63.367195  51.248614  "
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Seems to be a small bug in the homework\n",
    "# They said they kept the largest mean expression of duplicate probes, but\n",
    "# the homework code does not do this. For example, ZZZ3 has two probes\n",
    "# and the homework code keeps the smaller mean expression ILMN_2137536\n",
    "2 ** all_mixes[all_mixes['Gene symbol'] == 'ZZZ3'].drop('Gene symbol', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove duplicate rows with duplicate Gene symbols, keeping the row with the highest average expression\n",
    "# From https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4739640:\n",
    "# \"For non-Affymetrix platforms [this one is Illumina], genes mapping to >1 \n",
    "# probe were collapsed at the gene-level according to the probe with highest \n",
    "# mean expression across all samples.\"\n",
    "\n",
    "mean_expression_per_probe = all_mixes.drop('Gene symbol', axis=1, inplace=False).T.mean()\n",
    "best_duplicate_probe_index = pd.concat([mean_expression_per_probe, all_mixes['Gene symbol']], axis=1).sort_values(0).drop_duplicates(subset='Gene symbol', keep='last')\n",
    "deduped_mixes = all_mixes.loc[best_duplicate_probe_index.index].set_index('Gene symbol').sort_index()\n",
    "\n",
    "# Remove log transform to preserve linearity of mixing\n",
    "delogged = (2 ** deduped_mixes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [],
   "source": [
    "delogged.to_csv('./data/cell_mixes.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.11.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}