CCLE_T-Test.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "title: Drug Response\n",
    "path: CCLE/t-test\n",
    "authors:\n",
    "- kellrott\n",
    "- adamstruck\n",
    "tags:\n",
    "- ccle\n",
    "- drug response\n",
    "created_at: 2018-05-09\n",
    "updated_at: 2020-01-14\n",
    "tldr: Identify drug response seperation for CCLE cohorts defined by gene mutation\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import gripql\n",
    "import itertools\n",
    "import scipy.stats as stats\n",
    "\n",
    "conn = gripql.Connection(\"https://bmeg.io/api\", credential_file=\"bmeg_credentials.json\")\n",
    "G = conn.graph(\"rc5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find all of the samples in the [CTRP](https://portals.broadinstitute.org/ctrp/) Breast Cancer experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-01-14 13:36:52,216\t40 results received in 0 seconds\n"
     ]
    }
   ],
   "source": [
    "q = G.query().V(\"Program:CTRP\").out(\"projects\").out(\"cases\").\\\n",
    "    has(gripql.eq(\"cellline_attributes.Primary Disease\", \"Breast Cancer\")).distinct()\n",
    "all_cases = []\n",
    "for row in q:\n",
    "    all_cases.append(row.gid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the genes of interest, get Ensembl gene ids, from the HUGO symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "GENES = [\"PTEN\", \"TP53\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-01-14 13:36:56,003\t2 results received in 0 seconds\n"
     ]
    }
   ],
   "source": [
    "gene_ids = {}\n",
    "for i in G.query().V().hasLabel(\"Gene\").has(gripql.within(\"symbol\", GENES)):\n",
    "    gene_ids[i.data.symbol] = i.gid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CTRP doesn't have direct mutation calling, but rather they used the same cell lines as the CCLE, and we can use the mutation calls from that project. So use the `same_as` edge to identify the cases from CCLE that are the same as the ones in CTRP. Then use the mutations from those samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'TP53': 'ENSG00000141510', 'PTEN': 'ENSG00000171862'}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_ids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each of the genes, find the set of samples that have a mutation in that gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-01-14 13:39:36,629\t43 results received in 0 seconds\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ENSG00000141510 Positive Set: 29\n",
      "ENSG00000141510 Negative Set: 11\n",
      "ENSG00000171862 Positive Set: 9\n",
      "ENSG00000171862 Negative Set: 31\n"
     ]
    }
   ],
   "source": [
    "mut_cases = {}\n",
    "norm_cases = {}\n",
    "\n",
    "q = G.query().V(all_cases).as_(\"ctrp\").out(\"same_as\").has(gripql.eq(\"project_id\", \"Project:CCLE\"))\n",
    "q = q.out(\"samples\").out(\"aliquots\").out(\"somatic_callsets\")\n",
    "q = q.outE(\"alleles\").has(gripql.within(\"ensembl_gene\", list(gene_ids.values())))\n",
    "q = q.render({\"case\" : \"$ctrp._gid\", \"gene\" : \"$._data.ensembl_gene\"})\n",
    "\n",
    "for res in q:\n",
    "    mut_cases[res.gene] = mut_cases.get(res.gene, set()) | set([res.case])\n",
    "\n",
    "#get CCLE samples without mutation\n",
    "for i in gene_ids.values():\n",
    "    norm_cases[i] = list(set(all_cases).difference(mut_cases[i]))\n",
    "\n",
    "    print( \"%s Positive Set: %d\" % (i, len(mut_cases[i])) )\n",
    "    print( \"%s Negative Set: %d\" % (i, len(norm_cases[i])) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-01-14 13:39:55,140\t12,224 results received in 3 seconds\n",
      "[INFO]\t2020-01-14 13:39:56,594\t3,999 results received in 1 seconds\n"
     ]
    }
   ],
   "source": [
    "pos_response = {}\n",
    "for g in gene_ids.values():\n",
    "    pos_response[g] = {}\n",
    "    q = G.query().V(list(mut_cases[g])).as_(\"a\").out(\"samples\").out(\"aliquots\")\n",
    "    q = q.out(\"drug_response\").as_(\"a\").out(\"compounds\").as_(\"b\")\n",
    "    q = q.select([\"a\", \"b\"])    \n",
    "    for row in q:\n",
    "        v = row['a']['data']['aac']\n",
    "        compound = row['b']['gid']\n",
    "        if compound not in pos_response[g]:\n",
    "            pos_response[g][compound] = [ v ]\n",
    "        else:\n",
    "            pos_response[g][compound].append(v)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[INFO]\t2020-01-14 13:39:59,162\t5,025 results received in 1 seconds\n",
      "[INFO]\t2020-01-14 13:40:03,457\t13,250 results received in 4 seconds\n"
     ]
    }
   ],
   "source": [
    "neg_response = {}\n",
    "for g in gene_ids.values():\n",
    "    neg_response[g] = {}\n",
    "    q = G.query().V(list(norm_cases[g])).as_(\"a\").out(\"samples\").out(\"aliquots\")\n",
    "    q = q.out(\"drug_response\").as_(\"a\").out(\"compounds\").as_(\"b\")\n",
    "    q = q.select([\"a\", \"b\"])    \n",
    "    for row in q:\n",
    "        v = row['a']['data']['aac']\n",
    "        compound = row['b']['gid']\n",
    "        if compound not in neg_response[g]:\n",
    "            neg_response[g][compound] = [ v ]\n",
    "        else:\n",
    "            neg_response[g][compound].append(v)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "drugs = set(itertools.chain.from_iterable( i.keys() for i in pos_response.values() ))\n",
    "out = []\n",
    "for drug in drugs:\n",
    "    for g in gene_ids.values():\n",
    "        if drug in pos_response[g] and drug in neg_response[g]:\n",
    "            row = {\"drug\" : drug, \"mutation\" : g}\n",
    "            mut_values = pos_response[g][drug]\n",
    "            norm_values = neg_response[g][drug]\n",
    "            if len(mut_values) > 5 and len(norm_values) > 5:\n",
    "                s = stats.ttest_ind(mut_values, norm_values, equal_var=False)\n",
    "                row[\"t-statistic\"] = s.statistic\n",
    "                row[\"t-pvalue\"] = s.pvalue\n",
    "                s = stats.f_oneway(mut_values, norm_values)\n",
    "                row[\"a-statistic\"] = s.statistic\n",
    "                row[\"a-pvalue\"] = s.pvalue\n",
    "                out.append(row)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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>drug</th>\n",
       "      <th>mutation</th>\n",
       "      <th>t-statistic</th>\n",
       "      <th>t-pvalue</th>\n",
       "      <th>a-statistic</th>\n",
       "      <th>a-pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>277</td>\n",
       "      <td>Compound:CID216345</td>\n",
       "      <td>ENSG00000141510</td>\n",
       "      <td>-4.425335</td>\n",
       "      <td>0.000366</td>\n",
       "      <td>23.811104</td>\n",
       "      <td>0.000019</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>593</td>\n",
       "      <td>Compound:CID16736978</td>\n",
       "      <td>ENSG00000141510</td>\n",
       "      <td>-3.177435</td>\n",
       "      <td>0.006779</td>\n",
       "      <td>12.249536</td>\n",
       "      <td>0.001356</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>214</td>\n",
       "      <td>Compound:CID9967941</td>\n",
       "      <td>ENSG00000141510</td>\n",
       "      <td>-3.114875</td>\n",
       "      <td>0.006419</td>\n",
       "      <td>11.046079</td>\n",
       "      <td>0.002136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>31</td>\n",
       "      <td>Compound:CID9549297</td>\n",
       "      <td>ENSG00000141510</td>\n",
       "      <td>-2.573574</td>\n",
       "      <td>0.021992</td>\n",
       "      <td>10.077463</td>\n",
       "      <td>0.003071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>480</td>\n",
       "      <td>Compound:CID56949517</td>\n",
       "      <td>ENSG00000171862</td>\n",
       "      <td>2.235155</td>\n",
       "      <td>0.069114</td>\n",
       "      <td>10.480699</td>\n",
       "      <td>0.003638</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>808</td>\n",
       "      <td>Compound:CID5472285</td>\n",
       "      <td>ENSG00000141510</td>\n",
       "      <td>-2.711139</td>\n",
       "      <td>0.015176</td>\n",
       "      <td>9.369077</td>\n",
       "      <td>0.004037</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>566</td>\n",
       "      <td>Compound:NO_ONTOLOGY:CID9809715</td>\n",
       "      <td>ENSG00000141510</td>\n",
       "      <td>-2.492241</td>\n",
       "      <td>0.026354</td>\n",
       "      <td>8.810857</td>\n",
       "      <td>0.005631</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>237</td>\n",
       "      <td>Compound:CID451668</td>\n",
       "      <td>ENSG00000171862</td>\n",
       "      <td>2.030714</td>\n",
       "      <td>0.051533</td>\n",
       "      <td>7.157550</td>\n",
       "      <td>0.008627</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>95</td>\n",
       "      <td>Compound:CID54686904</td>\n",
       "      <td>ENSG00000171862</td>\n",
       "      <td>2.241852</td>\n",
       "      <td>0.067998</td>\n",
       "      <td>8.668157</td>\n",
       "      <td>0.009525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>656</td>\n",
       "      <td>Compound:CID11609586</td>\n",
       "      <td>ENSG00000171862</td>\n",
       "      <td>1.948182</td>\n",
       "      <td>0.066221</td>\n",
       "      <td>6.818910</td>\n",
       "      <td>0.011089</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                drug         mutation  t-statistic  t-pvalue  \\\n",
       "277               Compound:CID216345  ENSG00000141510    -4.425335  0.000366   \n",
       "593             Compound:CID16736978  ENSG00000141510    -3.177435  0.006779   \n",
       "214              Compound:CID9967941  ENSG00000141510    -3.114875  0.006419   \n",
       "31               Compound:CID9549297  ENSG00000141510    -2.573574  0.021992   \n",
       "480             Compound:CID56949517  ENSG00000171862     2.235155  0.069114   \n",
       "808              Compound:CID5472285  ENSG00000141510    -2.711139  0.015176   \n",
       "566  Compound:NO_ONTOLOGY:CID9809715  ENSG00000141510    -2.492241  0.026354   \n",
       "237               Compound:CID451668  ENSG00000171862     2.030714  0.051533   \n",
       "95              Compound:CID54686904  ENSG00000171862     2.241852  0.067998   \n",
       "656             Compound:CID11609586  ENSG00000171862     1.948182  0.066221   \n",
       "\n",
       "     a-statistic  a-pvalue  \n",
       "277    23.811104  0.000019  \n",
       "593    12.249536  0.001356  \n",
       "214    11.046079  0.002136  \n",
       "31     10.077463  0.003071  \n",
       "480    10.480699  0.003638  \n",
       "808     9.369077  0.004037  \n",
       "566     8.810857  0.005631  \n",
       "237     7.157550  0.008627  \n",
       "95      8.668157  0.009525  \n",
       "656     6.818910  0.011089  "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pandas.DataFrame(out, columns=[\"drug\", \"mutation\", \"t-statistic\", \"t-pvalue\", \"a-statistic\", \"a-pvalue\"]).sort_values(\"a-pvalue\").head(10)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}