3.build-hetnets/generate_msigdb_matrix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing MSigDB Gene Sets into Binary Matrix\n",
    "\n",
    "This notebook loads the full MSigDB gene set `.gmt` file (version 6.1) and outputs a binary, gene by gene set matrix indicating gene membership in the given gene set.\n",
    "\n",
    "**Note that we exclude gene sets with restrictive licences (KEGG, Biocarta, and The AAAS/STKE Cell Signaling Database)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import csv\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_template_matrix(msigdb_file, blacklist, checkblacklist=True):\n",
    "    \"\"\"\n",
    "    Retrieve all genes and pathways from given msigdb .gmt file\n",
    "    \n",
    "    Output:\n",
    "    sorted gene by pathways pandas dataframe. Entries indicate membership\n",
    "    \"\"\"\n",
    "    all_db_pathways = []\n",
    "    all_db_genes = []\n",
    "\n",
    "    # Get a set of all genes and all pathways in MSigDB (not blacklisted)\n",
    "    with open(msigdb_file, 'r') as msigdb_fh:\n",
    "        msigdb_reader = csv.reader(msigdb_fh, delimiter='\\t')\n",
    "\n",
    "        for row in msigdb_reader:\n",
    "            signature_name = row[0]\n",
    "            signature_genes = row[2:]\n",
    "            \n",
    "            if checkblacklist:\n",
    "                if signature_name.startswith(blacklist):\n",
    "                    continue\n",
    "\n",
    "            all_db_pathways.append(signature_name)\n",
    "            all_db_genes += signature_genes\n",
    "        \n",
    "    big_msigdb_df = pd.DataFrame(0, index=set(all_db_genes), columns=all_db_pathways)\n",
    "    big_msigdb_df = big_msigdb_df.sort_index()\n",
    "    big_msigdb_df = big_msigdb_df.T.sort_index().T\n",
    "    \n",
    "    # Loop through file again to populate dataframe. This is a fast implementation\n",
    "    with open(msigdb_file, 'r') as msigdb_fh:\n",
    "        msigdb_reader = csv.reader(msigdb_fh, delimiter='\\t')\n",
    "        for row in msigdb_reader:\n",
    "            signature_name = row[0]\n",
    "            signature_genes = row[2:]\n",
    "            if checkblacklist:\n",
    "                if signature_name.startswith(blacklist):\n",
    "                    continue\n",
    "\n",
    "            for gene in signature_genes:\n",
    "                big_msigdb_df.at[gene, signature_name] = 1\n",
    "\n",
    "    return big_msigdb_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store .gmt files\n",
    "full_msigdb_file = os.path.join('data', 'msigdb.v6.1.entrez.gmt')\n",
    "\n",
    "# Resources with restrictive licenses\n",
    "blacklist = ('KEGG', 'BIOCARTA', 'ST_')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process MSigDB gmt files into large matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(32431, 17356)\n",
      "CPU times: user 18 s, sys: 1.36 s, total: 19.4 s\n",
      "Wall time: 19.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "full_msigdb_df = make_template_matrix(full_msigdb_file, blacklist, checkblacklist=True)\n",
    "print(full_msigdb_df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10min 18s, sys: 52 ms, total: 10min 18s\n",
      "Wall time: 10min 16s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "full_msigdb_file = os.path.join('data', 'full_msigdb_binary_matrix.tsv.bz2')\n",
    "full_msigdb_df.to_csv(full_msigdb_file, sep='\\t', compression='bz2')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:interpret-compression]",
   "language": "python",
   "name": "conda-env-interpret-compression-py"
  },
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}