L1000_repurpoisng.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# L1000 repurpoing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as ss\n",
    "import h5py\n",
    "import csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input rank of zscore profile, and positions of up and down regulated genes to calculate ES score\n",
    "\n",
    "def calc_ES(zscorerank,qup,qdown):  \n",
    "    n = len(zscorerank)\n",
    "    nu = len(qup)\n",
    "    nd = len(qdown)\n",
    "\n",
    "    uptemp = np.array(sorted(zscorerank[qup]))\n",
    "    downtemp = np.array(sorted(zscorerank[qdown]))\n",
    "    upmax = max(np.array(range(1, nu+1))/nu - uptemp/n)\n",
    "    upmin = min(np.array(range(0, nu))/nu - (uptemp - 1)/n)\n",
    "    downmax = max(np.array(range(1, nd+1))/nd - downtemp/n)\n",
    "    downmin = min(np.array(range(0, nd))/nd - (downtemp - 1)/n)\n",
    "\n",
    "    upscore = (upmax if abs(upmax) > abs(upmin) else upmin)\n",
    "    downscore = (downmax if abs(downmax) > abs(downmin) else downmin)\n",
    "  \n",
    "    ES= ((downscore - upscore)/2 if(upscore*downscore < 0) else 0)\n",
    "    \n",
    "    return ES\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_csv(f):\n",
    "    l=[]\n",
    "    with open(f) as csv_file:\n",
    "        f_reader=csv.reader(csv_file)\n",
    "        for row in f_reader:\n",
    "            l.append(row)\n",
    "    return l\n",
    "\n",
    "def read_tsv(f):\n",
    "    l=[]\n",
    "    with open(f) as csv_file:\n",
    "        f_reader=csv.reader(csv_file,delimiter='\\t')\n",
    "        for row in f_reader:\n",
    "            l.append(row)\n",
    "    return l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import L1000 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path=\"/path/to/L1000_files/\"\n",
    "L1000_ph1 = h5py.File(data_path+'Bayesian_GSE92742_Level5_COMPZ_n361481x978.h5', 'r')\n",
    "L1000_ph2 = h5py.File(data_path+'Bayesian_GSE70138_Level5_COMPZ_n116218x978.h5', 'r')\n",
    "\n",
    "ph1_siginfo= read_tsv(data_path+\"GSE92742_Broad_LINCS_sig_info.txt\")\n",
    "ph2_siginfo= read_tsv(data_path+\"GSE70138_Broad_LINCS_sig_info.txt\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "ph1_sig_dict={}\n",
    "for i in range(len(ph1_siginfo)):\n",
    "    ph1_sig_dict[ph1_siginfo[i][0]]=ph1_siginfo[i]\n",
    "    \n",
    "ph2_sig_dict={}\n",
    "for i in range(len(ph2_siginfo)):\n",
    "    ph2_sig_dict[ph1_siginfo[i][0]]=ph2_siginfo[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import up/down regulated genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_path=\"/path/to/gene_files/\"\n",
    "up_genes= read_csv(gene_path+\"up_genes.csv\")[0]\n",
    "down_genes= read_csv(gene_path+\"down_genes.csv\")[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "ph1_sigs=[x.decode() for x in L1000_ph1['colid']]\n",
    "ph2_sigs=[x.decode() for x in L1000_ph2['colid']]\n",
    "lm_genes=[x.decode() for x in L1000_ph1['rowid']]\n",
    "\n",
    "ph1_data=np.array(L1000_ph1['data'])\n",
    "ph2_data=np.array(L1000_ph2['data'])\n",
    "\n",
    "lm_dict={}\n",
    "for i in range(len(lm_genes)):\n",
    "    lm_dict[lm_genes[i]]=i\n",
    "\n",
    "ph1_sig_dict={}\n",
    "for i in range(len(ph1_siginfo)):\n",
    "    ph1_sig_dict[ph1_siginfo[i][0]]=ph1_siginfo[i]\n",
    "    \n",
    "ph2_sig_dict={}\n",
    "for i in range(len(ph2_siginfo)):\n",
    "    ph2_sig_dict[ph2_siginfo[i][0]]=ph2_siginfo[i]\n",
    "\n",
    "\n",
    "L1000_up_genes=list(set([g.strip() for g in up_genes]) & set(lm_genes))\n",
    "L1000_down_genes=list(set([g.strip() for g in down_genes]) & set(lm_genes))\n",
    "\n",
    "qup=[lm_dict[g] for g in L1000_up_genes]\n",
    "qdown=[lm_dict[g] for g in L1000_up_genes]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate enrichment scrores for all L1000 profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "ph1_ranked=[ss.rankdata(i) for i in ph1_data]\n",
    "ph2_ranked=[ss.rankdata(i) for i in ph2_data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "ES_list1=[]\n",
    "for i in range(len(ph1_sigs)):\n",
    "    score= calc_ES(ph1_ranked[i],qup,qdown)\n",
    "    ES_list1.append(score)\n",
    "\n",
    "ES_list2=[]\n",
    "for i in range(len(ph2_sigs)):\n",
    "    score= calc_ES(ph2_ranked[i],qup,qdown)\n",
    "    ES_list2.append(score)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Top sigs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_N=10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "sort_index1 = np.argsort(ES_list1)\n",
    "sort_index2 = np.argsort(ES_list2)\n",
    "\n",
    "ph1_topsigs=[ph1_sigs[i] for i in sort_index1[-top_N:]]\n",
    "ph2_topsigs=[ph2_sigs[i] for i in sort_index2[-top_N:]]\n",
    "\n",
    "ph1_topinfo=[ph1_sig_dict[i] for i in ph1_topsigs]\n",
    "ph2_topinfo=[ph2_sig_dict[i] for i in ph2_topsigs]\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "7b4b4feff2f24a0f0a34464dbe537a36fda679851528fb8735cb41fa49dffb2d"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit ('base': conda)",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}