z_statistics_enrich.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1e1ab39b-f624-4d76-a987-d4825bdd7536",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0ab0d4da-0e86-4c0c-9ee9-c1f1b39ed8cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import pathlib\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95fb12e8-a334-4e1d-9796-169907154b82",
   "metadata": {},
   "source": [
    "### Classes and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e5ee936e-cc3b-4632-a62d-1bddb93cd22c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Synonym_Normal_Assumption():\n",
    "\n",
    "    def __init__(self, df_synonym, score_column=\"score\", se_column=\"SE\", one_sided=True):\n",
    "        \"\"\"\n",
    "        A class to perform a z-test on the scores of non-synonymous variants \n",
    "        compared to the distribution of all synonymous variants. This assumes that\n",
    "        standard errors for each score are known and uses a weighted mean\n",
    "        (by standard errors) to perform the z-Test.\n",
    "        \"\"\"\n",
    "        self.df_syn = df_synonym\n",
    "        self.__score = score_column\n",
    "        self.__se = se_column\n",
    "        self.__w = \"weights\"\n",
    "        self.one_sided = one_sided\n",
    "        self._fit()\n",
    "        \n",
    "    def _fit(self):\n",
    "        \"\"\"\n",
    "        calculate weights and means under the assumption that the synonyms constitute our\n",
    "        baseline normal distribution\n",
    "        \"\"\"\n",
    "        self.df_syn[self.__w] = 1/self.df_syn[self.__se]**2\n",
    "        self.weights_sum = self.df_syn[self.__w].sum()\n",
    "        self.weighted_mean = ((self.df_syn[self.__score]*self.df_syn[self.__w]).sum()) / self.weights_sum\n",
    "        self.sem = np.sqrt(1/((self.df_syn[self.__se]**2)**-1).sum())\n",
    "\n",
    "\n",
    "    def get_p_value_two_tailed_standard_normal(self, z_score):\n",
    "        \"\"\"calculate the p two-sided\"\"\"\n",
    "        return 2 * (1 - sp.stats.norm.cdf(abs(z_score)))\n",
    "    \n",
    "    def get_p_value_one_tailed_standard_normal(self, z_score):\n",
    "        \"\"\"calculate the p one-sided\"\"\"\n",
    "        return (1 - sp.stats.norm.cdf(abs(z_score)))\n",
    "\n",
    "    \n",
    "    def calc(self):\n",
    "        \"\"\"return a calc function that returns the raw p-value\"\"\"\n",
    "        def __calc(row):\n",
    "            z = self.calculate_z(row)\n",
    "            if self.one_sided:\n",
    "                return self.get_p_value_one_tailed_standard_normal(z)\n",
    "            else:\n",
    "                return self.get_p_value_two_tailed_standard_normal(z)\n",
    "        return __calc\n",
    "\n",
    "    def calculate_z(self, row):\n",
    "        \"\"\"Calculate z statistic\"\"\"\n",
    "        z = (row[self.__score] - self.weighted_mean) / self.calculate_sediff(row)\n",
    "        return z\n",
    "\n",
    "    def calculate_sediff(self, row):\n",
    "        \"\"\"Calculate z statistic\"\"\"\n",
    "        sediff = np.sqrt(self.sem**2 + row[self.__se]**2)\n",
    "        return sediff\n",
    "\n",
    "    def benjamini_hochberg_correction(self, p_values, name=\"adj. p-value(B.H.)\"):\n",
    "        padj = pd.Series(\n",
    "            sp.stats.false_discovery_control(p_values),\n",
    "            index=p_values.index,\n",
    "            name=name\n",
    "        )\n",
    "        return padj\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d09d1bd-fc4b-41f2-b622-18b39ae343b5",
   "metadata": {},
   "source": [
    "### Setttings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "501345de-1aca-476e-a86e-24ed61269006",
   "metadata": {},
   "outputs": [],
   "source": [
    "out_folder = Path(\"out\")\n",
    "input_table = out_folder / \"Exon5-8_Enrich2_RFS.tsv\"\n",
    "output_file = out_folder / \"Exon5-8_Enrich2_RFS_z.tsv\"\n",
    "\n",
    "id_column = \"mut_ID\"  # id of the variant sequence\n",
    "score_column=\"transformed_score_enrich2\"  # the score to test on \n",
    "se_column=\"transformed_SE_enrich2\"  # # the corresponding standard error\n",
    "syn_column = \"synonym\"  # column indicating if this is a synonymous mutation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caf01e00-a5d6-4224-862d-a6e521a57dde",
   "metadata": {},
   "source": [
    "### Z-Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "db042fe7-94bd-4d28-aceb-dcbc62a66eca",
   "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>hg38_genomic</th>\n",
       "      <th>hg38_cDNA</th>\n",
       "      <th>hg38_protein</th>\n",
       "      <th>read_count_dmso_1</th>\n",
       "      <th>read_count_dmso_2</th>\n",
       "      <th>read_count_dmso_3</th>\n",
       "      <th>read_count_n3a_1</th>\n",
       "      <th>read_count_n3a_2</th>\n",
       "      <th>read_count_n3a_3</th>\n",
       "      <th>read_count_donor</th>\n",
       "      <th>...</th>\n",
       "      <th>Rep2_score</th>\n",
       "      <th>Rep3_SE</th>\n",
       "      <th>Rep3_score</th>\n",
       "      <th>score_CIup</th>\n",
       "      <th>score_CIdown</th>\n",
       "      <th>transformed_score_enrich2</th>\n",
       "      <th>transformed_score_CIup</th>\n",
       "      <th>transformed_score_CIdown</th>\n",
       "      <th>transformed_SE_enrich2</th>\n",
       "      <th>synonym</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mut_ID</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",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NC_000017.11:g.7675248_7675250del</td>\n",
       "      <td>NM_000546.6:c.376-12_376-10del</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>842</td>\n",
       "      <td>614</td>\n",
       "      <td>624</td>\n",
       "      <td>8</td>\n",
       "      <td>58</td>\n",
       "      <td>237</td>\n",
       "      <td>3767</td>\n",
       "      <td>...</td>\n",
       "      <td>-2.758029</td>\n",
       "      <td>0.076236</td>\n",
       "      <td>-1.631491</td>\n",
       "      <td>-1.149392</td>\n",
       "      <td>-5.144778</td>\n",
       "      <td>-0.970377</td>\n",
       "      <td>0.298027</td>\n",
       "      <td>-2.238781</td>\n",
       "      <td>0.647145</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NC_000017.11:g.7675247_7675248del</td>\n",
       "      <td>NM_000546.6:c.376-10_376-9del</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>6506</td>\n",
       "      <td>2489</td>\n",
       "      <td>3260</td>\n",
       "      <td>1103</td>\n",
       "      <td>545</td>\n",
       "      <td>776</td>\n",
       "      <td>25047</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.924381</td>\n",
       "      <td>0.039932</td>\n",
       "      <td>-2.099547</td>\n",
       "      <td>-1.897810</td>\n",
       "      <td>-2.321020</td>\n",
       "      <td>-0.311525</td>\n",
       "      <td>-0.177169</td>\n",
       "      <td>-0.445880</td>\n",
       "      <td>0.068549</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NC_000017.11:g.7675249del</td>\n",
       "      <td>NM_000546.6:c.376-12del</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>4726</td>\n",
       "      <td>2135</td>\n",
       "      <td>2142</td>\n",
       "      <td>747</td>\n",
       "      <td>164</td>\n",
       "      <td>920</td>\n",
       "      <td>17449</td>\n",
       "      <td>...</td>\n",
       "      <td>-2.969792</td>\n",
       "      <td>0.039410</td>\n",
       "      <td>-1.509519</td>\n",
       "      <td>-1.450645</td>\n",
       "      <td>-3.109950</td>\n",
       "      <td>-0.420024</td>\n",
       "      <td>0.106751</td>\n",
       "      <td>-0.946799</td>\n",
       "      <td>0.268763</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NC_000017.11:g.7675248G&gt;A</td>\n",
       "      <td>NM_000546.6:c.376-12C&gt;T</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>4441</td>\n",
       "      <td>1576</td>\n",
       "      <td>2528</td>\n",
       "      <td>616</td>\n",
       "      <td>51</td>\n",
       "      <td>489</td>\n",
       "      <td>14982</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.827627</td>\n",
       "      <td>0.049381</td>\n",
       "      <td>-2.306705</td>\n",
       "      <td>-1.940106</td>\n",
       "      <td>-3.798452</td>\n",
       "      <td>-0.793988</td>\n",
       "      <td>-0.204025</td>\n",
       "      <td>-1.383952</td>\n",
       "      <td>0.301002</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>NC_000017.11:g.7675248G&gt;C</td>\n",
       "      <td>NM_000546.6:c.376-12C&gt;G</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>1680</td>\n",
       "      <td>937</td>\n",
       "      <td>622</td>\n",
       "      <td>375</td>\n",
       "      <td>2</td>\n",
       "      <td>113</td>\n",
       "      <td>5006</td>\n",
       "      <td>...</td>\n",
       "      <td>-6.333172</td>\n",
       "      <td>0.102064</td>\n",
       "      <td>-2.366648</td>\n",
       "      <td>-0.863307</td>\n",
       "      <td>-6.155259</td>\n",
       "      <td>-1.200349</td>\n",
       "      <td>0.479672</td>\n",
       "      <td>-2.880369</td>\n",
       "      <td>0.857153</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 55 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                             hg38_genomic                       hg38_cDNA  \\\n",
       "mut_ID                                                                      \n",
       "1       NC_000017.11:g.7675248_7675250del  NM_000546.6:c.376-12_376-10del   \n",
       "2       NC_000017.11:g.7675247_7675248del   NM_000546.6:c.376-10_376-9del   \n",
       "3               NC_000017.11:g.7675249del         NM_000546.6:c.376-12del   \n",
       "4               NC_000017.11:g.7675248G>A         NM_000546.6:c.376-12C>T   \n",
       "5               NC_000017.11:g.7675248G>C         NM_000546.6:c.376-12C>G   \n",
       "\n",
       "           hg38_protein  read_count_dmso_1  read_count_dmso_2  \\\n",
       "mut_ID                                                          \n",
       "1       NP_000537.3:p.?                842                614   \n",
       "2       NP_000537.3:p.?               6506               2489   \n",
       "3       NP_000537.3:p.?               4726               2135   \n",
       "4       NP_000537.3:p.?               4441               1576   \n",
       "5       NP_000537.3:p.?               1680                937   \n",
       "\n",
       "        read_count_dmso_3  read_count_n3a_1  read_count_n3a_2  \\\n",
       "mut_ID                                                          \n",
       "1                     624                 8                58   \n",
       "2                    3260              1103               545   \n",
       "3                    2142               747               164   \n",
       "4                    2528               616                51   \n",
       "5                     622               375                 2   \n",
       "\n",
       "        read_count_n3a_3  read_count_donor  ... Rep2_score   Rep3_SE  \\\n",
       "mut_ID                                      ...                        \n",
       "1                    237              3767  ...  -2.758029  0.076236   \n",
       "2                    776             25047  ...  -1.924381  0.039932   \n",
       "3                    920             17449  ...  -2.969792  0.039410   \n",
       "4                    489             14982  ...  -3.827627  0.049381   \n",
       "5                    113              5006  ...  -6.333172  0.102064   \n",
       "\n",
       "       Rep3_score score_CIup score_CIdown  transformed_score_enrich2  \\\n",
       "mut_ID                                                                 \n",
       "1       -1.631491  -1.149392    -5.144778                  -0.970377   \n",
       "2       -2.099547  -1.897810    -2.321020                  -0.311525   \n",
       "3       -1.509519  -1.450645    -3.109950                  -0.420024   \n",
       "4       -2.306705  -1.940106    -3.798452                  -0.793988   \n",
       "5       -2.366648  -0.863307    -6.155259                  -1.200349   \n",
       "\n",
       "       transformed_score_CIup  transformed_score_CIdown  \\\n",
       "mut_ID                                                    \n",
       "1                    0.298027                 -2.238781   \n",
       "2                   -0.177169                 -0.445880   \n",
       "3                    0.106751                 -0.946799   \n",
       "4                   -0.204025                 -1.383952   \n",
       "5                    0.479672                 -2.880369   \n",
       "\n",
       "       transformed_SE_enrich2 synonym  \n",
       "mut_ID                                 \n",
       "1                    0.647145   False  \n",
       "2                    0.068549   False  \n",
       "3                    0.268763   False  \n",
       "4                    0.301002   False  \n",
       "5                    0.857153   False  \n",
       "\n",
       "[5 rows x 55 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load an input dataframe with scores and standard errors \n",
    "\n",
    "df_enrich_rfs = pd.read_csv(input_table, sep=\"\\t\")\n",
    "df_enrich_rfs = df_enrich_rfs.set_index(id_column)\n",
    "df_enrich_rfs[\"synonym\"] = df_enrich_rfs[\"type_p\"] == \"syn\"\n",
    "display(df_enrich_rfs.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4113dcb9-175f-4f6c-a399-d5c19238e890",
   "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>mut_ID</th>\n",
       "      <th>hg38_genomic</th>\n",
       "      <th>hg38_cDNA</th>\n",
       "      <th>hg38_protein</th>\n",
       "      <th>read_count_dmso_1</th>\n",
       "      <th>read_count_dmso_2</th>\n",
       "      <th>read_count_dmso_3</th>\n",
       "      <th>read_count_n3a_1</th>\n",
       "      <th>read_count_n3a_2</th>\n",
       "      <th>read_count_n3a_3</th>\n",
       "      <th>...</th>\n",
       "      <th>score_CIup</th>\n",
       "      <th>score_CIdown</th>\n",
       "      <th>transformed_score_enrich2</th>\n",
       "      <th>transformed_score_CIup</th>\n",
       "      <th>transformed_score_CIdown</th>\n",
       "      <th>transformed_SE_enrich2</th>\n",
       "      <th>synonym</th>\n",
       "      <th>weights</th>\n",
       "      <th>p (one-sided, all)</th>\n",
       "      <th>p_adjusted (one-sided, all)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>377</th>\n",
       "      <td>1</td>\n",
       "      <td>NC_000017.11:g.7675248_7675250del</td>\n",
       "      <td>NM_000546.6:c.376-12_376-10del</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>842</td>\n",
       "      <td>614</td>\n",
       "      <td>624</td>\n",
       "      <td>8</td>\n",
       "      <td>58</td>\n",
       "      <td>237</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.149392</td>\n",
       "      <td>-5.144778</td>\n",
       "      <td>-0.970377</td>\n",
       "      <td>0.298027</td>\n",
       "      <td>-2.238781</td>\n",
       "      <td>0.647145</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.492765</td>\n",
       "      <td>0.493714</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>378</th>\n",
       "      <td>2</td>\n",
       "      <td>NC_000017.11:g.7675247_7675248del</td>\n",
       "      <td>NM_000546.6:c.376-10_376-9del</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>6506</td>\n",
       "      <td>2489</td>\n",
       "      <td>3260</td>\n",
       "      <td>1103</td>\n",
       "      <td>545</td>\n",
       "      <td>776</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.897810</td>\n",
       "      <td>-2.321020</td>\n",
       "      <td>-0.311525</td>\n",
       "      <td>-0.177169</td>\n",
       "      <td>-0.445880</td>\n",
       "      <td>0.068549</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>379</th>\n",
       "      <td>3</td>\n",
       "      <td>NC_000017.11:g.7675249del</td>\n",
       "      <td>NM_000546.6:c.376-12del</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>4726</td>\n",
       "      <td>2135</td>\n",
       "      <td>2142</td>\n",
       "      <td>747</td>\n",
       "      <td>164</td>\n",
       "      <td>920</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.450645</td>\n",
       "      <td>-3.109950</td>\n",
       "      <td>-0.420024</td>\n",
       "      <td>0.106751</td>\n",
       "      <td>-0.946799</td>\n",
       "      <td>0.268763</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.022544</td>\n",
       "      <td>0.027076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>380</th>\n",
       "      <td>4</td>\n",
       "      <td>NC_000017.11:g.7675248G&gt;A</td>\n",
       "      <td>NM_000546.6:c.376-12C&gt;T</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>4441</td>\n",
       "      <td>1576</td>\n",
       "      <td>2528</td>\n",
       "      <td>616</td>\n",
       "      <td>51</td>\n",
       "      <td>489</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.940106</td>\n",
       "      <td>-3.798452</td>\n",
       "      <td>-0.793988</td>\n",
       "      <td>-0.204025</td>\n",
       "      <td>-1.383952</td>\n",
       "      <td>0.301002</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.292202</td>\n",
       "      <td>0.309629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>381</th>\n",
       "      <td>5</td>\n",
       "      <td>NC_000017.11:g.7675248G&gt;C</td>\n",
       "      <td>NM_000546.6:c.376-12C&gt;G</td>\n",
       "      <td>NP_000537.3:p.?</td>\n",
       "      <td>1680</td>\n",
       "      <td>937</td>\n",
       "      <td>622</td>\n",
       "      <td>375</td>\n",
       "      <td>2</td>\n",
       "      <td>113</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.863307</td>\n",
       "      <td>-6.155259</td>\n",
       "      <td>-1.200349</td>\n",
       "      <td>0.479672</td>\n",
       "      <td>-2.880369</td>\n",
       "      <td>0.857153</td>\n",
       "      <td>False</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.388977</td>\n",
       "      <td>0.400334</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 59 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     mut_ID                       hg38_genomic  \\\n",
       "377       1  NC_000017.11:g.7675248_7675250del   \n",
       "378       2  NC_000017.11:g.7675247_7675248del   \n",
       "379       3          NC_000017.11:g.7675249del   \n",
       "380       4          NC_000017.11:g.7675248G>A   \n",
       "381       5          NC_000017.11:g.7675248G>C   \n",
       "\n",
       "                          hg38_cDNA     hg38_protein  read_count_dmso_1  \\\n",
       "377  NM_000546.6:c.376-12_376-10del  NP_000537.3:p.?                842   \n",
       "378   NM_000546.6:c.376-10_376-9del  NP_000537.3:p.?               6506   \n",
       "379         NM_000546.6:c.376-12del  NP_000537.3:p.?               4726   \n",
       "380         NM_000546.6:c.376-12C>T  NP_000537.3:p.?               4441   \n",
       "381         NM_000546.6:c.376-12C>G  NP_000537.3:p.?               1680   \n",
       "\n",
       "     read_count_dmso_2  read_count_dmso_3  read_count_n3a_1  read_count_n3a_2  \\\n",
       "377                614                624                 8                58   \n",
       "378               2489               3260              1103               545   \n",
       "379               2135               2142               747               164   \n",
       "380               1576               2528               616                51   \n",
       "381                937                622               375                 2   \n",
       "\n",
       "     read_count_n3a_3  ...  score_CIup score_CIdown transformed_score_enrich2  \\\n",
       "377               237  ...   -1.149392    -5.144778                 -0.970377   \n",
       "378               776  ...   -1.897810    -2.321020                 -0.311525   \n",
       "379               920  ...   -1.450645    -3.109950                 -0.420024   \n",
       "380               489  ...   -1.940106    -3.798452                 -0.793988   \n",
       "381               113  ...   -0.863307    -6.155259                 -1.200349   \n",
       "\n",
       "    transformed_score_CIup transformed_score_CIdown transformed_SE_enrich2  \\\n",
       "377               0.298027                -2.238781               0.647145   \n",
       "378              -0.177169                -0.445880               0.068549   \n",
       "379               0.106751                -0.946799               0.268763   \n",
       "380              -0.204025                -1.383952               0.301002   \n",
       "381               0.479672                -2.880369               0.857153   \n",
       "\n",
       "     synonym weights  p (one-sided, all) p_adjusted (one-sided, all)  \n",
       "377    False     NaN            0.492765                    0.493714  \n",
       "378    False     NaN                 0.0                         0.0  \n",
       "379    False     NaN            0.022544                    0.027076  \n",
       "380    False     NaN            0.292202                    0.309629  \n",
       "381    False     NaN            0.388977                    0.400334  \n",
       "\n",
       "[5 rows x 59 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# divide the variants into synonymous and non-synonymous mutations\n",
    "concat = []\n",
    "df_synonym = df_enrich_rfs.loc[df_enrich_rfs[syn_column], :].copy()\n",
    "df_non_synonym = df_enrich_rfs.loc[(~df_enrich_rfs[syn_column]), :].copy()\n",
    "\n",
    "# Perform z-test\n",
    "syn = Synonym_Normal_Assumption(\n",
    "        df_synonym, \n",
    "        score_column=score_column, \n",
    "        se_column=se_column\n",
    "    )\n",
    "df_non_synonym[\"p (one-sided, all)\"] = df_non_synonym.apply(syn.calc(), axis=1)\n",
    "\n",
    "# calcualte Benjamini-Hochberg correction\n",
    "df_non_synonym[\"p_adjusted (one-sided, all)\"] = syn.benjamini_hochberg_correction(\n",
    "    df_non_synonym[\"p (one-sided, all)\"]\n",
    ")\n",
    "\n",
    "# concat the synonymous and non-synoymous variants dataframes\n",
    "df_synonym[\"p (one-sided, all)\"] = \"NA\"\n",
    "df_synonym[\"p_adjusted (one-sided, all)\"] = \"NA\"\n",
    "concat.append(df_synonym)\n",
    "concat.append(df_non_synonym)\n",
    "df_enrich_rfs_z = pd.concat(concat)\n",
    "\n",
    "# reset index and sort\n",
    "df_enrich_rfs_z = df_enrich_rfs_z.reset_index()\n",
    "df_enrich_rfs_z = df_enrich_rfs_z.sort_values(id_column)\n",
    "\n",
    "# display result frame\n",
    "display(df_enrich_rfs_z.head())\n",
    "\n",
    "# write to file\n",
    "df_enrich_rfs_z.to_csv(output_file, index=False, sep=\"\\t\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}