GEO_data_retrieval.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "15c877f6-d17d-4f37-bac0-055a415c3ba0",
   "metadata": {
    "tags": []
   },
   "source": [
    "This notebook provides the functionality to retrieve gene expression data from the Gene Expression Omnibus (GEO), apply data normalization techniques, and generate an output table in the .tsv format."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfa44387-8c80-41c7-859e-1e57f37d7671",
   "metadata": {},
   "source": [
    "Example will be shown on Ulloa-Montoya GSE35640, GPL570 sample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8ed2d69-704f-4be9-96cd-edaaf6b2b8db",
   "metadata": {},
   "source": [
    "# Import of python base packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3b5fca7a-9113-45a0-bbcb-cdb1db08178f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%matplotlib inline\n",
    "%config IPCompleter.use_jedi = False\n",
    "\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import pathlib\n",
    "import subprocess\n",
    "import logging\n",
    "import csv\n",
    "import matplotlib.pyplot as plt\n",
    "import umap\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import SVG\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "%config InlineBackend.figure_format = 'png'\n",
    "plt.rcParams['pdf.fonttype'] = 'truetype'\n",
    "plt.rcParams['svg.fonttype'] = 'none'\n",
    "plt.rcParams['figure.dpi'] = 120"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fad0b441-09a1-49af-98e0-35a3936b2dad",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%matplotlib inline\n",
    "%config IPCompleter.use_jedi = False\n",
    "\n",
    "import pathlib\n",
    "import subprocess\n",
    "import logging\n",
    "import os\n",
    "import pandas as pd\n",
    "import csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e46bab67-0f3c-402f-87a3-2d3db738c900",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from portraits.mapping import get_gs_for_probes_from_3col,get_expressions_for_gs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "eaa2346c-7421-4aa3-9b37-c15faf93b756",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8928db8d-fdce-4b11-87ca-57187ab24b90",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a89260e-74f2-4395-b5df-9de3a5fb3be7",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Get the example data from GEO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51b5d784-1f07-49e6-8cdb-84157af90040",
   "metadata": {},
   "source": [
    "If you want to use another sample please change the GSE and PLATFORM variables to your desired samples GSE and PLATFORM values in GEO correspondingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "81c220ac-d0fc-4d03-8699-4e12dec04c3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "GSE = 'GSE35640'\n",
    "PLATFORM = 'GPL570'\n",
    "EXPRESSION_MATRIX = 'Cohorts/TCGA_input_expressions/expression_from_GEO.tsv'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c803755-fe38-47c8-b438-b4fa37fab871",
   "metadata": {},
   "source": [
    "Initialize and create a temporary directory for the CEL files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8681c081-1416-47b2-a6eb-8f7aad18ef73",
   "metadata": {},
   "outputs": [],
   "source": [
    "current_dir = pathlib.Path().parent.absolute()\n",
    "dir_to_process = str(current_dir / 'TMPDIR')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e9697d13-840a-4122-a45b-60daee0508ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(dir_to_process):\n",
    "     os.mkdir(dir_to_process)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7dc6c27e-ff19-4d8a-84b6-5ba89ff99bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.devnull, \"w\") as f:\n",
    "    subprocess.run([\n",
    "        'wget',\n",
    "        f'ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE35nnn/{GSE}/suppl/{GSE}_RAW.tar'\n",
    "    ],stdout=f, stderr=subprocess.STDOUT );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ba3c61ee-d311-496d-8ce6-203bcba2018e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tar: Skipping to next header\n",
      "tar: Exiting with failure status due to previous errors\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['rm', 'GSE35640_RAW.tar'], returncode=0)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subprocess.run([\n",
    "    'tar',\n",
    "    '-xf',\n",
    "    f'{GSE}_RAW.tar',\n",
    "    '-C', dir_to_process\n",
    "])\n",
    "subprocess.run([\n",
    "    'rm',\n",
    "    f'{GSE}_RAW.tar'\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1f619c48-a5b0-4d11-a4ba-5f166519661f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['GSM872328_MAGE008_sample_1.CEL.gz']"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.listdir(dir_to_process)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8c78800-88e0-43c4-8ea9-9178ae00cd46",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Extracting expression values from CEL file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c486edb3-a7b9-43a8-950a-6bdda6822f47",
   "metadata": {},
   "source": [
    "For affy arrays without special probes, use RMA\n",
    "For GPL570/GPL96, use gcrma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1efa13ab-8169-40fa-9fce-d404f699d59a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-02T13:25:21.735530Z",
     "start_time": "2018-08-02T13:25:16.965609Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: BiocGenerics\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘BiocGenerics’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    IQR, mad, sd, var, xtabs\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:base’:\n",
      "\n",
      "    anyDuplicated, append, as.data.frame, basename, cbind, colnames,\n",
      "    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,\n",
      "    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,\n",
      "    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,\n",
      "    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,\n",
      "    union, unique, unsplit, which.max, which.min\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: Biobase\n",
      "\n",
      "R[write to console]: Welcome to Bioconductor\n",
      "\n",
      "    Vignettes contain introductory material; view with\n",
      "    'browseVignettes()'. To cite Bioconductor, see\n",
      "    'citation(\"Biobase\")', and for packages 'citation(\"pkgname\")'.\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: AnnotationDbi\n",
      "\n",
      "R[write to console]: Loading required package: stats4\n",
      "\n",
      "R[write to console]: Loading required package: IRanges\n",
      "\n",
      "R[write to console]: Loading required package: S4Vectors\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘S4Vectors’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:base’:\n",
      "\n",
      "    expand.grid, I, unname\n",
      "\n",
      "\n",
      "R[write to console]: Loading required package: XML\n",
      "\n",
      "R[write to console]: \n",
      "Attaching package: ‘XML’\n",
      "\n",
      "\n",
      "R[write to console]: The following object is masked from ‘package:tools’:\n",
      "\n",
      "    toHTML\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "# Load required R packages\n",
    "library(affy)\n",
    "library(annotate)\n",
    "library(gcrma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a9f47bd1-1757-44a8-afa4-61590b1b5c7e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-02T14:25:38.950721Z",
     "start_time": "2018-08-02T14:19:31.888815Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adjusting for optical effect.................................................................Done.\n",
      "Computing affinities.Done.\n",
      "Adjusting for non-specific binding.................................................................Done.\n",
      "Normalizing\n",
      "Calculating Expression\n"
     ]
    }
   ],
   "source": [
    "%%R -i dir_to_process -o normalized_expression\n",
    "\n",
    "# Bulk read cell files\n",
    "raw_expression <- ReadAffy(celfile.path = dir_to_process)\n",
    "\n",
    "# Normalize expression using RMA\n",
    "rma_normalized <- gcrma(raw_expression)\n",
    "\n",
    "# Retrieve expressions from dataset\n",
    "normalized_expression <- as.data.frame(exprs(rma_normalized))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d422bec2-59b9-4ce8-b860-677cd2eb4057",
   "metadata": {
    "scrolled": true
   },
   "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>GSM872328_MAGE008_sample_1.CEL.gz</th>\n",
       "      <th>GSM872329_MAGE008_sample_2.CEL.gz</th>\n",
       "      <th>GSM872330_MAGE008_sample_3.CEL.gz</th>\n",
       "      <th>GSM872331_MAGE008_sample_4.CEL.gz</th>\n",
       "      <th>GSM872332_MAGE008_sample_5.CEL.gz</th>\n",
       "      <th>GSM872333_MAGE008_sample_6.CEL.gz</th>\n",
       "      <th>GSM872334_MAGE008_sample_7.CEL.gz</th>\n",
       "      <th>GSM872335_MAGE008_sample_8.CEL.gz</th>\n",
       "      <th>GSM872336_MAGE008_sample_9.CEL.gz</th>\n",
       "      <th>GSM872337_MAGE008_sample_10.CEL.gz</th>\n",
       "      <th>...</th>\n",
       "      <th>GSM872383_MAGE008_sample_56.CEL.gz</th>\n",
       "      <th>GSM872384_MAGE008_sample_57.CEL.gz</th>\n",
       "      <th>GSM872385_MAGE008_sample_58.CEL.gz</th>\n",
       "      <th>GSM872386_MAGE008_sample_59.CEL.gz</th>\n",
       "      <th>GSM872387_MAGE008_sample_60.CEL.gz</th>\n",
       "      <th>GSM872388_MAGE008_sample_61.CEL.gz</th>\n",
       "      <th>GSM872389_MAGE008_sample_62.CEL.gz</th>\n",
       "      <th>GSM872390_MAGE008_sample_63.CEL.gz</th>\n",
       "      <th>GSM872391_MAGE008_sample_64.CEL.gz</th>\n",
       "      <th>GSM872392_MAGE008_sample_65.CEL.gz</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1007_s_at</th>\n",
       "      <td>9.816090</td>\n",
       "      <td>6.474644</td>\n",
       "      <td>9.113994</td>\n",
       "      <td>7.044742</td>\n",
       "      <td>8.719946</td>\n",
       "      <td>4.852191</td>\n",
       "      <td>7.190953</td>\n",
       "      <td>8.636192</td>\n",
       "      <td>5.623183</td>\n",
       "      <td>7.433451</td>\n",
       "      <td>...</td>\n",
       "      <td>6.287400</td>\n",
       "      <td>7.017124</td>\n",
       "      <td>7.533525</td>\n",
       "      <td>7.042195</td>\n",
       "      <td>8.809050</td>\n",
       "      <td>8.198514</td>\n",
       "      <td>8.457781</td>\n",
       "      <td>7.530108</td>\n",
       "      <td>4.724734</td>\n",
       "      <td>9.222943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1053_at</th>\n",
       "      <td>5.664499</td>\n",
       "      <td>6.540040</td>\n",
       "      <td>5.569006</td>\n",
       "      <td>5.745150</td>\n",
       "      <td>4.690936</td>\n",
       "      <td>5.923728</td>\n",
       "      <td>5.700713</td>\n",
       "      <td>5.120923</td>\n",
       "      <td>5.394368</td>\n",
       "      <td>6.704893</td>\n",
       "      <td>...</td>\n",
       "      <td>5.992869</td>\n",
       "      <td>6.754348</td>\n",
       "      <td>5.943736</td>\n",
       "      <td>6.333925</td>\n",
       "      <td>6.063615</td>\n",
       "      <td>6.360003</td>\n",
       "      <td>6.073347</td>\n",
       "      <td>7.141974</td>\n",
       "      <td>7.191548</td>\n",
       "      <td>5.611221</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117_at</th>\n",
       "      <td>8.394219</td>\n",
       "      <td>6.083483</td>\n",
       "      <td>6.587381</td>\n",
       "      <td>6.092671</td>\n",
       "      <td>3.411626</td>\n",
       "      <td>3.791668</td>\n",
       "      <td>5.667734</td>\n",
       "      <td>4.965991</td>\n",
       "      <td>4.475886</td>\n",
       "      <td>4.443469</td>\n",
       "      <td>...</td>\n",
       "      <td>3.296701</td>\n",
       "      <td>5.638913</td>\n",
       "      <td>3.205921</td>\n",
       "      <td>3.471334</td>\n",
       "      <td>3.682211</td>\n",
       "      <td>8.971476</td>\n",
       "      <td>4.343425</td>\n",
       "      <td>5.738228</td>\n",
       "      <td>4.807307</td>\n",
       "      <td>10.331267</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>121_at</th>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.345001</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.388071</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.454035</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>3.128954</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>...</td>\n",
       "      <td>2.320004</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.411101</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.380696</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.882576</td>\n",
       "      <td>2.434246</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1255_g_at</th>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>...</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.348901</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 65 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           GSM872328_MAGE008_sample_1.CEL.gz  \\\n",
       "1007_s_at                           9.816090   \n",
       "1053_at                             5.664499   \n",
       "117_at                              8.394219   \n",
       "121_at                              2.390202   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872329_MAGE008_sample_2.CEL.gz  \\\n",
       "1007_s_at                           6.474644   \n",
       "1053_at                             6.540040   \n",
       "117_at                              6.083483   \n",
       "121_at                              2.345001   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872330_MAGE008_sample_3.CEL.gz  \\\n",
       "1007_s_at                           9.113994   \n",
       "1053_at                             5.569006   \n",
       "117_at                              6.587381   \n",
       "121_at                              2.390202   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872331_MAGE008_sample_4.CEL.gz  \\\n",
       "1007_s_at                           7.044742   \n",
       "1053_at                             5.745150   \n",
       "117_at                              6.092671   \n",
       "121_at                              2.388071   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872332_MAGE008_sample_5.CEL.gz  \\\n",
       "1007_s_at                           8.719946   \n",
       "1053_at                             4.690936   \n",
       "117_at                              3.411626   \n",
       "121_at                              2.390202   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872333_MAGE008_sample_6.CEL.gz  \\\n",
       "1007_s_at                           4.852191   \n",
       "1053_at                             5.923728   \n",
       "117_at                              3.791668   \n",
       "121_at                              2.390202   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872334_MAGE008_sample_7.CEL.gz  \\\n",
       "1007_s_at                           7.190953   \n",
       "1053_at                             5.700713   \n",
       "117_at                              5.667734   \n",
       "121_at                              2.454035   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872335_MAGE008_sample_8.CEL.gz  \\\n",
       "1007_s_at                           8.636192   \n",
       "1053_at                             5.120923   \n",
       "117_at                              4.965991   \n",
       "121_at                              2.390202   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872336_MAGE008_sample_9.CEL.gz  \\\n",
       "1007_s_at                           5.623183   \n",
       "1053_at                             5.394368   \n",
       "117_at                              4.475886   \n",
       "121_at                              3.128954   \n",
       "1255_g_at                           2.223272   \n",
       "\n",
       "           GSM872337_MAGE008_sample_10.CEL.gz  ...  \\\n",
       "1007_s_at                            7.433451  ...   \n",
       "1053_at                              6.704893  ...   \n",
       "117_at                               4.443469  ...   \n",
       "121_at                               2.390202  ...   \n",
       "1255_g_at                            2.223272  ...   \n",
       "\n",
       "           GSM872383_MAGE008_sample_56.CEL.gz  \\\n",
       "1007_s_at                            6.287400   \n",
       "1053_at                              5.992869   \n",
       "117_at                               3.296701   \n",
       "121_at                               2.320004   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872384_MAGE008_sample_57.CEL.gz  \\\n",
       "1007_s_at                            7.017124   \n",
       "1053_at                              6.754348   \n",
       "117_at                               5.638913   \n",
       "121_at                               2.390202   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872385_MAGE008_sample_58.CEL.gz  \\\n",
       "1007_s_at                            7.533525   \n",
       "1053_at                              5.943736   \n",
       "117_at                               3.205921   \n",
       "121_at                               2.411101   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872386_MAGE008_sample_59.CEL.gz  \\\n",
       "1007_s_at                            7.042195   \n",
       "1053_at                              6.333925   \n",
       "117_at                               3.471334   \n",
       "121_at                               2.390202   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872387_MAGE008_sample_60.CEL.gz  \\\n",
       "1007_s_at                            8.809050   \n",
       "1053_at                              6.063615   \n",
       "117_at                               3.682211   \n",
       "121_at                               2.390202   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872388_MAGE008_sample_61.CEL.gz  \\\n",
       "1007_s_at                            8.198514   \n",
       "1053_at                              6.360003   \n",
       "117_at                               8.971476   \n",
       "121_at                               2.390202   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872389_MAGE008_sample_62.CEL.gz  \\\n",
       "1007_s_at                            8.457781   \n",
       "1053_at                              6.073347   \n",
       "117_at                               4.343425   \n",
       "121_at                               2.380696   \n",
       "1255_g_at                            2.348901   \n",
       "\n",
       "           GSM872390_MAGE008_sample_63.CEL.gz  \\\n",
       "1007_s_at                            7.530108   \n",
       "1053_at                              7.141974   \n",
       "117_at                               5.738228   \n",
       "121_at                               2.390202   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872391_MAGE008_sample_64.CEL.gz  \\\n",
       "1007_s_at                            4.724734   \n",
       "1053_at                              7.191548   \n",
       "117_at                               4.807307   \n",
       "121_at                               2.882576   \n",
       "1255_g_at                            2.223272   \n",
       "\n",
       "           GSM872392_MAGE008_sample_65.CEL.gz  \n",
       "1007_s_at                            9.222943  \n",
       "1053_at                              5.611221  \n",
       "117_at                              10.331267  \n",
       "121_at                               2.434246  \n",
       "1255_g_at                            2.223272  \n",
       "\n",
       "[5 rows x 65 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normalized_expression.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "50c6cb3c-ebb9-4208-9795-f217d76ba3ef",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-02T14:25:38.956494Z",
     "start_time": "2018-08-02T14:25:38.952824Z"
    },
    "tags": []
   },
   "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>GSM872328</th>\n",
       "      <th>GSM872329</th>\n",
       "      <th>GSM872330</th>\n",
       "      <th>GSM872331</th>\n",
       "      <th>GSM872332</th>\n",
       "      <th>GSM872333</th>\n",
       "      <th>GSM872334</th>\n",
       "      <th>GSM872335</th>\n",
       "      <th>GSM872336</th>\n",
       "      <th>GSM872337</th>\n",
       "      <th>...</th>\n",
       "      <th>GSM872383</th>\n",
       "      <th>GSM872384</th>\n",
       "      <th>GSM872385</th>\n",
       "      <th>GSM872386</th>\n",
       "      <th>GSM872387</th>\n",
       "      <th>GSM872388</th>\n",
       "      <th>GSM872389</th>\n",
       "      <th>GSM872390</th>\n",
       "      <th>GSM872391</th>\n",
       "      <th>GSM872392</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1007_s_at</th>\n",
       "      <td>9.816090</td>\n",
       "      <td>6.474644</td>\n",
       "      <td>9.113994</td>\n",
       "      <td>7.044742</td>\n",
       "      <td>8.719946</td>\n",
       "      <td>4.852191</td>\n",
       "      <td>7.190953</td>\n",
       "      <td>8.636192</td>\n",
       "      <td>5.623183</td>\n",
       "      <td>7.433451</td>\n",
       "      <td>...</td>\n",
       "      <td>6.287400</td>\n",
       "      <td>7.017124</td>\n",
       "      <td>7.533525</td>\n",
       "      <td>7.042195</td>\n",
       "      <td>8.809050</td>\n",
       "      <td>8.198514</td>\n",
       "      <td>8.457781</td>\n",
       "      <td>7.530108</td>\n",
       "      <td>4.724734</td>\n",
       "      <td>9.222943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1053_at</th>\n",
       "      <td>5.664499</td>\n",
       "      <td>6.540040</td>\n",
       "      <td>5.569006</td>\n",
       "      <td>5.745150</td>\n",
       "      <td>4.690936</td>\n",
       "      <td>5.923728</td>\n",
       "      <td>5.700713</td>\n",
       "      <td>5.120923</td>\n",
       "      <td>5.394368</td>\n",
       "      <td>6.704893</td>\n",
       "      <td>...</td>\n",
       "      <td>5.992869</td>\n",
       "      <td>6.754348</td>\n",
       "      <td>5.943736</td>\n",
       "      <td>6.333925</td>\n",
       "      <td>6.063615</td>\n",
       "      <td>6.360003</td>\n",
       "      <td>6.073347</td>\n",
       "      <td>7.141974</td>\n",
       "      <td>7.191548</td>\n",
       "      <td>5.611221</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>117_at</th>\n",
       "      <td>8.394219</td>\n",
       "      <td>6.083483</td>\n",
       "      <td>6.587381</td>\n",
       "      <td>6.092671</td>\n",
       "      <td>3.411626</td>\n",
       "      <td>3.791668</td>\n",
       "      <td>5.667734</td>\n",
       "      <td>4.965991</td>\n",
       "      <td>4.475886</td>\n",
       "      <td>4.443469</td>\n",
       "      <td>...</td>\n",
       "      <td>3.296701</td>\n",
       "      <td>5.638913</td>\n",
       "      <td>3.205921</td>\n",
       "      <td>3.471334</td>\n",
       "      <td>3.682211</td>\n",
       "      <td>8.971476</td>\n",
       "      <td>4.343425</td>\n",
       "      <td>5.738228</td>\n",
       "      <td>4.807307</td>\n",
       "      <td>10.331267</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>121_at</th>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.345001</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.388071</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.454035</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>3.128954</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>...</td>\n",
       "      <td>2.320004</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.411101</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.380696</td>\n",
       "      <td>2.390202</td>\n",
       "      <td>2.882576</td>\n",
       "      <td>2.434246</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1255_g_at</th>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>...</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.348901</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "      <td>2.223272</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 65 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           GSM872328  GSM872329  GSM872330  GSM872331  GSM872332  GSM872333  \\\n",
       "1007_s_at   9.816090   6.474644   9.113994   7.044742   8.719946   4.852191   \n",
       "1053_at     5.664499   6.540040   5.569006   5.745150   4.690936   5.923728   \n",
       "117_at      8.394219   6.083483   6.587381   6.092671   3.411626   3.791668   \n",
       "121_at      2.390202   2.345001   2.390202   2.388071   2.390202   2.390202   \n",
       "1255_g_at   2.223272   2.223272   2.223272   2.223272   2.223272   2.223272   \n",
       "\n",
       "           GSM872334  GSM872335  GSM872336  GSM872337  ...  GSM872383  \\\n",
       "1007_s_at   7.190953   8.636192   5.623183   7.433451  ...   6.287400   \n",
       "1053_at     5.700713   5.120923   5.394368   6.704893  ...   5.992869   \n",
       "117_at      5.667734   4.965991   4.475886   4.443469  ...   3.296701   \n",
       "121_at      2.454035   2.390202   3.128954   2.390202  ...   2.320004   \n",
       "1255_g_at   2.223272   2.223272   2.223272   2.223272  ...   2.223272   \n",
       "\n",
       "           GSM872384  GSM872385  GSM872386  GSM872387  GSM872388  GSM872389  \\\n",
       "1007_s_at   7.017124   7.533525   7.042195   8.809050   8.198514   8.457781   \n",
       "1053_at     6.754348   5.943736   6.333925   6.063615   6.360003   6.073347   \n",
       "117_at      5.638913   3.205921   3.471334   3.682211   8.971476   4.343425   \n",
       "121_at      2.390202   2.411101   2.390202   2.390202   2.390202   2.380696   \n",
       "1255_g_at   2.223272   2.223272   2.223272   2.223272   2.223272   2.348901   \n",
       "\n",
       "           GSM872390  GSM872391  GSM872392  \n",
       "1007_s_at   7.530108   4.724734   9.222943  \n",
       "1053_at     7.141974   7.191548   5.611221  \n",
       "117_at      5.738228   4.807307  10.331267  \n",
       "121_at      2.390202   2.882576   2.434246  \n",
       "1255_g_at   2.223272   2.223272   2.223272  \n",
       "\n",
       "[5 rows x 65 columns]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Trim names to make the table more readable.\n",
    "\n",
    "normalized_expression.columns = normalized_expression.columns.to_series().apply(lambda x: x.split('_')[0]).values\n",
    "normalized_expression.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cf9c95de-3bd6-4f39-acf0-67400dff402c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['rm', '-r', 'TMPDIR/'], returncode=0)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Delete unnecessary files\n",
    "\n",
    "subprocess.run(['rm',  '-r','TMPDIR/'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70302308-2e72-44de-82f1-bd0401c17458",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Converting probe ids to HUGO gene symbols"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23919ffd-b3b0-426e-b64a-8d0d78b190ae",
   "metadata": {
    "tags": []
   },
   "source": [
    "To download the SOFT file manually, follow these steps (we use the Ulloya-Montoya sample as an example):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4287ff5-c398-46d6-af22-0e8a4e9ca0f1",
   "metadata": {},
   "source": [
    "* Go to https://www.ncbi.nlm.nih.gov/geo/\n",
    "* Type the GLP platform number – GLP570 – in the search bar\n",
    "* Go to https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570\n",
    "* Find the Annotation SOFT table at the bottom of the page\n",
    "* Click the button to download the SOFT annotation table for your sample (GPL570.annot.gz in our example)\n",
    "* Upload the file to the environment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08a80c9c-61a8-4930-954d-ac4bc1216768",
   "metadata": {},
   "source": [
    "To download SOFT annotation table from jupyter notebook, change the PLATFORM value to appropriate GPL platform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e69c6890-cde0-4297-be33-3eb10aacceea",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(os.devnull, \"w\") as f:\n",
    "    subprocess.run([\n",
    "        'wget',\n",
    "        f'ftp://ftp.ncbi.nlm.nih.gov/geo/platforms/GPLnnn/{PLATFORM}/annot/{PLATFORM}.annot.gz'\n",
    "    ], stdout=f, stderr=subprocess.STDOUT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c7dad83-c3ca-4c64-8ef9-2fd312ff8127",
   "metadata": {},
   "source": [
    "Once you’ve downloaded the SOFT file, extract the data from it. Unzip the SOFT file to get GPL570.annot.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c3e66f8e-18fb-49a0-9f0e-b834f07faaaa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['gunzip', 'GPL570.annot.gz'], returncode=0)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subprocess.run(['gunzip', f'{PLATFORM}.annot.gz'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12c1c8ec-d7e7-46f3-9409-7c1dc94ffa6e",
   "metadata": {},
   "source": [
    "Remove the header from the SOFT file to avoid any problems in further processes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "64e7aa93-ae67-4d2e-aa2b-7b0e43ddafc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=\"sed '1,/^ID/d' GPL570.annot > GPL570.beheaded.annot\", returncode=0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subprocess.run(\"sed '1,/^ID/d' GPL570.annot > GPL570.beheaded.annot\", shell = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5b62391-8b55-4bcc-9e43-b70c631b210e",
   "metadata": {},
   "source": [
    "Turn the SOFT file into a 3-column file.<br>\n",
    "\n",
    "1st column: probe id<br>\n",
    "2nd column: gene symbol column (as is with '///')<br>\n",
    "3rd column: entrez id (not needed for the study)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a62b7d78-13a2-4c9a-96da-235a78a8a17e",
   "metadata": {},
   "source": [
    "Subsetting SOFT file to have 3 columns "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "2525161e-22cc-416f-9eb8-86ec40b839d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_SOFT_annotations = pd.read_csv(f'{PLATFORM}.beheaded.annot', sep = '\\t', header = None)\n",
    "gene_SOFT_annotations = gene_SOFT_annotations.iloc[:, [0, 2, 3]]\n",
    "\n",
    "# Rename columns\n",
    "gene_SOFT_annotations = gene_SOFT_annotations.rename(columns = {2: 1, 3 : 2})\n",
    "gene_SOFT_annotations.to_csv(f'{PLATFORM}.3col', sep = '\\t',  index=False, quoting=csv.QUOTE_NONNUMERIC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "05c4e4ee-5687-403d-a0ea-d4c0a715f15a",
   "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>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1007_s_at</td>\n",
       "      <td>MIR4640///DDR1</td>\n",
       "      <td>100616237///780</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1053_at</td>\n",
       "      <td>RFC2</td>\n",
       "      <td>5982</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>117_at</td>\n",
       "      <td>HSPA6</td>\n",
       "      <td>3310</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>121_at</td>\n",
       "      <td>PAX8</td>\n",
       "      <td>7849</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1255_g_at</td>\n",
       "      <td>GUCA1A</td>\n",
       "      <td>2978</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           0               1                2\n",
       "0  1007_s_at  MIR4640///DDR1  100616237///780\n",
       "1    1053_at            RFC2             5982\n",
       "2     117_at           HSPA6             3310\n",
       "3     121_at            PAX8             7849\n",
       "4  1255_g_at          GUCA1A             2978"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_SOFT_annotations.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01313fa8-9287-4ae5-b718-0c4178528b9b",
   "metadata": {},
   "source": [
    "Delete all unnecessary files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "9f85fdc9-ef2d-4301-a51f-297df2ce762b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['rm', 'GPL570.beheaded.annot'], returncode=0)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "subprocess.run(['rm', f'{PLATFORM}.annot'])\n",
    "subprocess.run(['rm', f'{PLATFORM}.beheaded.annot'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6a65c011-7a90-465e-89ac-a5439289c713",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-02T14:25:50.137216Z",
     "start_time": "2018-08-02T14:25:49.716242Z"
    }
   },
   "outputs": [],
   "source": [
    "probes_gs_dict = get_gs_for_probes_from_3col(f'{PLATFORM}.3col', normalized_expression.index.tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "bf12dade-60f2-4132-a4b9-c624a8e01036",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1007_s_at    [MIR4640, DDR1]\n",
       "1053_at               [RFC2]\n",
       "117_at               [HSPA6]\n",
       "121_at                [PAX8]\n",
       "1255_g_at           [GUCA1A]\n",
       "1294_at      [MIR5193, UBA7]\n",
       "1316_at               [THRA]\n",
       "1320_at             [PTPN21]\n",
       "1405_i_at             [CCL5]\n",
       "1431_at             [CYP2E1]\n",
       "dtype: object"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(probes_gs_dict).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4f01fdfd-e8ff-4efb-8c20-102e8385e7f0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-02T14:26:46.823094Z",
     "start_time": "2018-08-02T14:25:52.167425Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "series = pd.Series(probes_gs_dict)\n",
    "annotated_expression = get_expressions_for_gs(series, normalized_expression, 'max').T.sort_index()\n",
    "\n",
    "annotated_expression.to_csv(EXPRESSION_MATRIX, sep='\\t', index=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "TME_env",
   "language": "python",
   "name": "tme_env"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}