variant_annotation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8c86d39f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np \n",
    "from scipy.stats import pearsonr\n",
    "from sklearn.metrics import classification_report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c45b841e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm, fisher_exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "523cc558",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K562_cov_atac.tsv\n",
      "All data, Pearson R2 0.09699208465104389\n",
      "FDR < 0.05, data, Pearson R2 0.2612582119502353\n",
      "Data without thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref               145   72\n",
      "alt               101  154\n",
      "SignificanceResult(statistic=3.0706820682068208, pvalue=4.654851873535244e-09)\n",
      "FDR < 0.05, model pval, data, Pearson R2 0.33434609309767915\n",
      "Data with thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref                96   38\n",
      "alt                47   93\n",
      "SignificanceResult(statistic=4.998880179171333, pvalue=2.3731087609957486e-10)\n",
      "K562_cov_dnase.tsv\n",
      "All data, Pearson R2 0.16271177747975757\n",
      "FDR < 0.05, data, Pearson R2 0.3487605002191143\n",
      "Data without thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref               786  399\n",
      "alt               435  809\n",
      "SignificanceResult(statistic=3.663607294097312, pvalue=1.4892206224325915e-54)\n",
      "FDR < 0.05, model pval, data, Pearson R2 0.42887225533406903\n",
      "Data with thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref               540  201\n",
      "alt               227  539\n",
      "SignificanceResult(statistic=6.3791176277204285, pvalue=5.245912765448729e-65)\n",
      "HepG2_cov_atac.tsv\n",
      "All data, Pearson R2 0.1533985991633921\n",
      "FDR < 0.05, data, Pearson R2 0.4500990870666304\n",
      "Data without thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref               137   63\n",
      "alt                40  102\n",
      "SignificanceResult(statistic=5.545238095238095, pvalue=1.890105468756956e-13)\n",
      "FDR < 0.05, model pval, data, Pearson R2 0.582462314647176\n",
      "Data with thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref                86   22\n",
      "alt                11   73\n",
      "SignificanceResult(statistic=25.94214876033058, pvalue=4.084574764117459e-21)\n",
      "HepG2_cov_dnase.tsv\n",
      "All data, Pearson R2 0.22395264375147161\n",
      "FDR < 0.05, data, Pearson R2 0.4266021682211808\n",
      "Data without thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref                 9    4\n",
      "alt                 6   14\n",
      "SignificanceResult(statistic=5.25, pvalue=0.03770534087794813)\n",
      "FDR < 0.05, model pval, data, Pearson R2 0.5812475038442471\n",
      "Data with thresholding model predictions\n",
      "pred_pref_allele  ref  alt\n",
      "pref_allele               \n",
      "ref                 6    1\n",
      "alt                 2   10\n",
      "SignificanceResult(statistic=30.0, pvalue=0.006271334444709059)\n"
     ]
    }
   ],
   "source": [
    "for path in [\"K562_cov_atac.tsv\", \"K562_cov_dnase.tsv\", \"HepG2_cov_atac.tsv\", \"HepG2_cov_dnase.tsv\"]:\n",
    "    print(path)\n",
    "    data = pd.read_csv(f\"../human_legnet/{path}\", sep=\"\\t\")\n",
    "    data['ref_pred'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].mean(axis=1)\n",
    "    data['alt_pred'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].mean(axis=1)\n",
    "    N =  data.columns.str.startswith(\"pred_ref\").sum()\n",
    "\n",
    "    var_ref = data.loc[:, data.columns.str.startswith(\"pred_ref\")].var(axis=1).mean()\n",
    "    var_alt = data.loc[:, data.columns.str.startswith(\"pred_alt\")].var(axis=1).mean()\n",
    "    var = var_ref + var_alt\n",
    "    std = np.sqrt(var / N)\n",
    "    #se = std * 1.64485\n",
    "\n",
    "    print(\"All data, Pearson R2\", pearsonr(data['score'], data['ref_pred'] - data['alt_pred'])[0])\n",
    "    \n",
    "    data['pred_score'] = data['ref_pred'] - data['alt_pred']\n",
    "    \n",
    "    data['pred_pval'] = norm.cdf(-abs(data['pred_score']) / std)\n",
    "    pred_significant =  data['pred_pval'] < 0.05\n",
    "    \n",
    "    \n",
    "    data['pred_pref_allele'] = (data['ref_pred'] - data['alt_pred']).apply(lambda x: \"ref\" if x > 0 else \"alt\")\n",
    "    \n",
    "    mask = data['fdr_comb_pval'] < 0.05\n",
    "\n",
    "    print(\"FDR < 0.05, data, Pearson R2\", pearsonr(data[mask]['score'], data[mask]['ref_pred'] - data[mask]['alt_pred'])[0])\n",
    "    \n",
    "    data_sig = data[mask].copy()\n",
    "    crosstab = pd.crosstab(index=data_sig['pref_allele'], columns=data_sig['pred_pref_allele'])\n",
    "    crosstab = crosstab.loc[['ref', 'alt'], ['ref', 'alt']]\n",
    "    print(\"Data without thresholding model predictions\")\n",
    "    print(crosstab)\n",
    "    print(fisher_exact(crosstab))\n",
    "    \n",
    "    \n",
    "    mask = np.logical_and(data['fdr_comb_pval'] < 0.05, pred_significant)\n",
    "    data_sig = data[mask].copy()\n",
    "    print(\"FDR < 0.05, model pval, data, Pearson R2\", pearsonr(data_sig['score'], \n",
    "                                                               data_sig['ref_pred'] - data_sig['alt_pred'])[0])\n",
    "    \n",
    "    \n",
    "    crosstab = pd.crosstab(index=data_sig['pref_allele'], columns=data_sig['pred_pref_allele'])\n",
    "    crosstab = crosstab.loc[['ref', 'alt'], ['ref', 'alt']]\n",
    "    print(\"Data with thresholding model predictions\")\n",
    "    print(crosstab)\n",
    "    print(fisher_exact(crosstab))\n",
    "    \n",
    "    final_data =  data[['chr',\n",
    "      'start',\n",
    "      'end', \n",
    "      'ref',\n",
    "      'alt', \n",
    "      'ref_pred',\n",
    "      'alt_pred',\n",
    "      'fdr_comb_pval',\n",
    "      'pred_pval',\n",
    "      'score',\n",
    "      'pred_score',\n",
    "      'pref_allele',\n",
    "      'pred_pref_allele']]\n",
    "    \n",
    "    name = path.replace(\".tsv\", \"_analysis.tsv\")\n",
    "    #final_data.to_csv(name, sep=\"\\t\", index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "id": "d7888237",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "274"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.logical_and(final_data['fdr_comb_pval'] < 0.05, final_data['pred_pval'] < 0.05).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "id": "e6851e57",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "189"
      ]
     },
     "execution_count": 207,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " 96 + 93"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "52c96e28",
   "metadata": {},
   "outputs": [],
   "source": [
    "eps=1e-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "aac046ef",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K562_asb.tsv\n",
      "All data, Pearson R2 0.07088269406962822\n",
      "FDR < 0.05, data, Pearson R2 0.16694896021607136\n",
      "Data without thresholding model predictions\n",
      "pred_pref_allele   ref   alt\n",
      "pref_allele                 \n",
      "ref               7212  5560\n",
      "alt               6590  8106\n",
      "SignificanceResult(statistic=1.5955194812283708, pvalue=1.9783535510627323e-82)\n",
      "FDR < 0.05, model pval, data, Pearson R2 0.21331963355095973\n",
      "Data with thresholding model predictions\n",
      "pred_pref_allele   ref   alt\n",
      "pref_allele                 \n",
      "ref               4290  2868\n",
      "alt               3311  4628\n",
      "SignificanceResult(statistic=2.0907991492792504, pvalue=2.2678725750157555e-111)\n",
      "HepG2_asb.tsv\n",
      "All data, Pearson R2 0.08018506978382701\n",
      "FDR < 0.05, data, Pearson R2 0.26180842634035567\n",
      "Data without thresholding model predictions\n",
      "pred_pref_allele   ref   alt\n",
      "pref_allele                 \n",
      "ref               2959  2086\n",
      "alt               1868  2667\n",
      "SignificanceResult(statistic=2.025241438282339, pvalue=1.7651682662523761e-65)\n",
      "FDR < 0.05, model pval, data, Pearson R2 0.3468084770651573\n",
      "Data with thresholding model predictions\n",
      "pred_pref_allele   ref   alt\n",
      "pref_allele                 \n",
      "ref               1661   887\n",
      "alt                781  1518\n",
      "SignificanceResult(statistic=3.639709735300189, pvalue=4.941428772403959e-106)\n"
     ]
    }
   ],
   "source": [
    "for path in [\"K562_asb.tsv\", \"HepG2_asb.tsv\"]:\n",
    "    print(path)\n",
    "    data = pd.read_csv(f\"../human_legnet/{path}\", sep=\"\\t\")\n",
    "    data['ref_pred'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].mean(axis=1)\n",
    "    data['alt_pred'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].mean(axis=1)\n",
    "    data['pref_allele'] = data[['fdrp_bh_ref',\n",
    "           'fdrp_bh_alt']].apply(lambda x: 'ref' if x.fdrp_bh_ref < x.fdrp_bh_alt else \"alt\", axis=1)\n",
    "    data['score'] = np.log2(data['min_fdr'] + eps) * (data['pref_allele'].apply(lambda x: -1 if x=='ref' else 1))\n",
    "\n",
    "\n",
    "    N =  data.columns.str.startswith(\"pred_ref\").sum()\n",
    "\n",
    "    var_ref = data.loc[:, data.columns.str.startswith(\"pred_ref\")].var(axis=1).mean()\n",
    "    var_alt = data.loc[:, data.columns.str.startswith(\"pred_alt\")].var(axis=1).mean()\n",
    "    var = var_ref + var_alt\n",
    "    std = np.sqrt(var / N)\n",
    "    #se = std * 1.64485\n",
    "\n",
    "    print(\"All data, Pearson R2\", pearsonr(data['score'], data['ref_pred'] - data['alt_pred'])[0])\n",
    "    \n",
    "    data['pred_score'] = data['ref_pred'] - data['alt_pred']\n",
    "    \n",
    "    data['pred_pval'] = norm.cdf(-abs(data['pred_score']) / std)\n",
    "    pred_significant =  data['pred_pval'] < 0.05\n",
    "    \n",
    "    \n",
    "    data['pred_pref_allele'] = (data['ref_pred'] - data['alt_pred']).apply(lambda x: \"ref\" if x > 0 else \"alt\")\n",
    "    \n",
    "    mask = data['min_fdr'] < 0.05\n",
    "\n",
    "    print(\"FDR < 0.05, data, Pearson R2\", pearsonr(data[mask]['score'], data[mask]['ref_pred'] - data[mask]['alt_pred'])[0])\n",
    "    \n",
    "    data_sig = data[mask].copy()\n",
    "    crosstab = pd.crosstab(index=data_sig['pref_allele'], columns=data_sig['pred_pref_allele'])\n",
    "    crosstab = crosstab.loc[['ref', 'alt'], ['ref', 'alt']]\n",
    "    print(\"Data without thresholding model predictions\")\n",
    "    print(crosstab)\n",
    "    print(fisher_exact(crosstab))\n",
    "    \n",
    "    data['fdr_comb_pval'] = data['min_fdr'] \n",
    "    data['start'] = data['pos']\n",
    "    data['end'] = data['pos'] + 1\n",
    "    \n",
    "    mask = np.logical_and(data['fdr_comb_pval'] < 0.05, pred_significant)\n",
    "    data_sig = data[mask].copy()\n",
    "    print(\"FDR < 0.05, model pval, data, Pearson R2\", pearsonr(data_sig['score'], \n",
    "                                                               data_sig['ref_pred'] - data_sig['alt_pred'])[0])\n",
    "    \n",
    "    \n",
    "    crosstab = pd.crosstab(index=data_sig['pref_allele'], columns=data_sig['pred_pref_allele'])\n",
    "    crosstab = crosstab.loc[['ref', 'alt'], ['ref', 'alt']]\n",
    "    print(\"Data with thresholding model predictions\")\n",
    "    print(crosstab)\n",
    "    print(fisher_exact(crosstab))\n",
    "    \n",
    "    final_data =  data[['chr',\n",
    "      'start',\n",
    "      'end', \n",
    "      'ref',\n",
    "      'alt', \n",
    "      'ref_pred',\n",
    "      'alt_pred',\n",
    "      'fdr_comb_pval',\n",
    "      'pred_pval',\n",
    "      'score',\n",
    "      'pred_score',\n",
    "      'pref_allele',\n",
    "      'pred_pref_allele']]\n",
    "    \n",
    "    name = path.replace(\".tsv\", \"_analysis.tsv\")\n",
    "    final_data.to_csv(name, sep=\"\\t\", index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "id": "fb635e7b",
   "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>chr</th>\n",
       "      <th>pos</th>\n",
       "      <th>ID</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>repeat_type</th>\n",
       "      <th>n_peak_calls</th>\n",
       "      <th>n_peak_callers</th>\n",
       "      <th>mean_BAD</th>\n",
       "      <th>mean_SNP_per_segment</th>\n",
       "      <th>...</th>\n",
       "      <th>pred_ref_test6_val3_0_rev</th>\n",
       "      <th>pred_alt_test6_val3_0_forw</th>\n",
       "      <th>pred_alt_test6_val3_0_rev</th>\n",
       "      <th>ref_pred</th>\n",
       "      <th>alt_pred</th>\n",
       "      <th>pref_allele</th>\n",
       "      <th>pred_score</th>\n",
       "      <th>pred_pval</th>\n",
       "      <th>pred_pref_allele</th>\n",
       "      <th>fdr_comb_pval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>770502</td>\n",
       "      <td>rs72631875</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>SINE</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>1105.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.40330</td>\n",
       "      <td>-0.4248</td>\n",
       "      <td>-0.4020</td>\n",
       "      <td>-0.407058</td>\n",
       "      <td>-0.421958</td>\n",
       "      <td>no</td>\n",
       "      <td>0.014899</td>\n",
       "      <td>1.582573e-02</td>\n",
       "      <td>ref</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>888410</td>\n",
       "      <td>rs143626389</td>\n",
       "      <td>G</td>\n",
       "      <td>A</td>\n",
       "      <td>Low_complexity</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.10</td>\n",
       "      <td>805.9</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.21510</td>\n",
       "      <td>-0.3235</td>\n",
       "      <td>-0.1935</td>\n",
       "      <td>-0.263535</td>\n",
       "      <td>-0.243948</td>\n",
       "      <td>alt</td>\n",
       "      <td>-0.019587</td>\n",
       "      <td>2.365418e-03</td>\n",
       "      <td>alt</td>\n",
       "      <td>2.364817e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>903551</td>\n",
       "      <td>rs7523690</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>LTR</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.03</td>\n",
       "      <td>1075.7</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.52600</td>\n",
       "      <td>-0.4888</td>\n",
       "      <td>-0.5566</td>\n",
       "      <td>-0.508631</td>\n",
       "      <td>-0.528689</td>\n",
       "      <td>no</td>\n",
       "      <td>0.020058</td>\n",
       "      <td>1.909637e-03</td>\n",
       "      <td>ref</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>905373</td>\n",
       "      <td>rs4970382</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>LINE</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.09</td>\n",
       "      <td>1305.8</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.22730</td>\n",
       "      <td>-0.2622</td>\n",
       "      <td>-0.2125</td>\n",
       "      <td>-0.113151</td>\n",
       "      <td>-0.012569</td>\n",
       "      <td>no</td>\n",
       "      <td>-0.100583</td>\n",
       "      <td>5.555454e-48</td>\n",
       "      <td>alt</td>\n",
       "      <td>8.868125e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>908025</td>\n",
       "      <td>rs11516185</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.17</td>\n",
       "      <td>1579.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.56300</td>\n",
       "      <td>-0.5110</td>\n",
       "      <td>-0.5547</td>\n",
       "      <td>-0.557515</td>\n",
       "      <td>-0.559646</td>\n",
       "      <td>no</td>\n",
       "      <td>0.002131</td>\n",
       "      <td>3.793201e-01</td>\n",
       "      <td>no</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>327934</th>\n",
       "      <td>chrX</td>\n",
       "      <td>155179766</td>\n",
       "      <td>rs6567793</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>LINE</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.50</td>\n",
       "      <td>235.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.26340</td>\n",
       "      <td>-0.2803</td>\n",
       "      <td>-0.2600</td>\n",
       "      <td>-0.261345</td>\n",
       "      <td>-0.254826</td>\n",
       "      <td>alt</td>\n",
       "      <td>-0.006519</td>\n",
       "      <td>1.735502e-01</td>\n",
       "      <td>no</td>\n",
       "      <td>5.104239e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>327935</th>\n",
       "      <td>chrX</td>\n",
       "      <td>155613005</td>\n",
       "      <td>rs89278</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.50</td>\n",
       "      <td>166.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.31180</td>\n",
       "      <td>0.5140</td>\n",
       "      <td>0.4045</td>\n",
       "      <td>0.489073</td>\n",
       "      <td>0.537469</td>\n",
       "      <td>no</td>\n",
       "      <td>-0.048396</td>\n",
       "      <td>1.480480e-12</td>\n",
       "      <td>alt</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>327936</th>\n",
       "      <td>chrX</td>\n",
       "      <td>155662569</td>\n",
       "      <td>rs553678</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1.50</td>\n",
       "      <td>235.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.20910</td>\n",
       "      <td>-0.2487</td>\n",
       "      <td>-0.2081</td>\n",
       "      <td>-0.212276</td>\n",
       "      <td>-0.203760</td>\n",
       "      <td>no</td>\n",
       "      <td>-0.008516</td>\n",
       "      <td>1.096976e-01</td>\n",
       "      <td>no</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>327937</th>\n",
       "      <td>chrX</td>\n",
       "      <td>155881414</td>\n",
       "      <td>rs28413172</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.25</td>\n",
       "      <td>547.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.11664</td>\n",
       "      <td>0.0511</td>\n",
       "      <td>0.1456</td>\n",
       "      <td>-0.007651</td>\n",
       "      <td>0.033403</td>\n",
       "      <td>no</td>\n",
       "      <td>-0.041055</td>\n",
       "      <td>1.601459e-09</td>\n",
       "      <td>alt</td>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>327938</th>\n",
       "      <td>chrX</td>\n",
       "      <td>156009119</td>\n",
       "      <td>rs111066114</td>\n",
       "      <td>G</td>\n",
       "      <td>C</td>\n",
       "      <td>Simple_repeat</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.00</td>\n",
       "      <td>928.0</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.19520</td>\n",
       "      <td>-0.2210</td>\n",
       "      <td>-0.1765</td>\n",
       "      <td>-0.047072</td>\n",
       "      <td>-0.040781</td>\n",
       "      <td>no</td>\n",
       "      <td>-0.006292</td>\n",
       "      <td>1.821055e-01</td>\n",
       "      <td>no</td>\n",
       "      <td>9.143509e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>327939 rows × 386 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         chr        pos           ID ref alt     repeat_type  n_peak_calls  \\\n",
       "0       chr1     770502   rs72631875   G   A            SINE             0   \n",
       "1       chr1     888410  rs143626389   G   A  Low_complexity             0   \n",
       "2       chr1     903551    rs7523690   A   C             LTR             0   \n",
       "3       chr1     905373    rs4970382   T   C            LINE             0   \n",
       "4       chr1     908025   rs11516185   A   G             NaN             0   \n",
       "...      ...        ...          ...  ..  ..             ...           ...   \n",
       "327934  chrX  155179766    rs6567793   A   G            LINE             0   \n",
       "327935  chrX  155613005      rs89278   C   T             NaN             0   \n",
       "327936  chrX  155662569     rs553678   A   G             NaN             0   \n",
       "327937  chrX  155881414   rs28413172   C   T             NaN             0   \n",
       "327938  chrX  156009119  rs111066114   G   C   Simple_repeat             0   \n",
       "\n",
       "        n_peak_callers  mean_BAD  mean_SNP_per_segment  ...  \\\n",
       "0                    0      1.00                1105.0  ...   \n",
       "1                    0      1.10                 805.9  ...   \n",
       "2                    0      1.03                1075.7  ...   \n",
       "3                    0      1.09                1305.8  ...   \n",
       "4                    0      1.17                1579.0  ...   \n",
       "...                ...       ...                   ...  ...   \n",
       "327934               0      1.50                 235.0  ...   \n",
       "327935               0      2.50                 166.0  ...   \n",
       "327936               0      1.50                 235.0  ...   \n",
       "327937               0      2.25                 547.0  ...   \n",
       "327938               0      2.00                 928.0  ...   \n",
       "\n",
       "        pred_ref_test6_val3_0_rev  pred_alt_test6_val3_0_forw  \\\n",
       "0                        -0.40330                     -0.4248   \n",
       "1                        -0.21510                     -0.3235   \n",
       "2                        -0.52600                     -0.4888   \n",
       "3                        -0.22730                     -0.2622   \n",
       "4                        -0.56300                     -0.5110   \n",
       "...                           ...                         ...   \n",
       "327934                   -0.26340                     -0.2803   \n",
       "327935                    0.31180                      0.5140   \n",
       "327936                   -0.20910                     -0.2487   \n",
       "327937                    0.11664                      0.0511   \n",
       "327938                   -0.19520                     -0.2210   \n",
       "\n",
       "        pred_alt_test6_val3_0_rev  ref_pred  alt_pred  pref_allele  \\\n",
       "0                         -0.4020 -0.407058 -0.421958           no   \n",
       "1                         -0.1935 -0.263535 -0.243948          alt   \n",
       "2                         -0.5566 -0.508631 -0.528689           no   \n",
       "3                         -0.2125 -0.113151 -0.012569           no   \n",
       "4                         -0.5547 -0.557515 -0.559646           no   \n",
       "...                           ...       ...       ...          ...   \n",
       "327934                    -0.2600 -0.261345 -0.254826          alt   \n",
       "327935                     0.4045  0.489073  0.537469           no   \n",
       "327936                    -0.2081 -0.212276 -0.203760           no   \n",
       "327937                     0.1456 -0.007651  0.033403           no   \n",
       "327938                    -0.1765 -0.047072 -0.040781           no   \n",
       "\n",
       "        pred_score     pred_pval  pred_pref_allele  fdr_comb_pval  \n",
       "0         0.014899  1.582573e-02               ref   1.000000e+00  \n",
       "1        -0.019587  2.365418e-03               alt   2.364817e-03  \n",
       "2         0.020058  1.909637e-03               ref   1.000000e+00  \n",
       "3        -0.100583  5.555454e-48               alt   8.868125e-02  \n",
       "4         0.002131  3.793201e-01                no   1.000000e+00  \n",
       "...            ...           ...               ...            ...  \n",
       "327934   -0.006519  1.735502e-01                no   5.104239e-07  \n",
       "327935   -0.048396  1.480480e-12               alt   1.000000e+00  \n",
       "327936   -0.008516  1.096976e-01                no   1.000000e+00  \n",
       "327937   -0.041055  1.601459e-09               alt   1.000000e+00  \n",
       "327938   -0.006292  1.821055e-01                no   9.143509e-01  \n",
       "\n",
       "[327939 rows x 386 columns]"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.2678e-111"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50d11f77",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "id": "dd4edca2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All data, Pearson R2 0.07088269406962822\n",
      "FDR < 0.05, data, Pearson R2 0.16694896021607136\n",
      "pred_pref      ref    alt      no\n",
      "pref_allele                      \n",
      "ref           4290   2868    5614\n",
      "alt           3311   4628    6757\n",
      "no           62173  61472  176826\n",
      "All data, Pearson R2 0.08018506978382701\n",
      "FDR < 0.05, data, Pearson R2 0.26180842634035567\n",
      "pred_pref      ref    alt      no\n",
      "pref_allele                      \n",
      "ref           1661    887    2497\n",
      "alt            781   1518    2236\n",
      "no           51233  50806  174677\n"
     ]
    }
   ],
   "source": [
    "for path in [\"K562_asb.tsv\", \"HepG2_asb.tsv\"]:\n",
    "\n",
    "    data = pd.read_csv(f\"../human_legnet/{path}\", sep=\"\\t\")\n",
    "    data['ref_pred'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].mean(axis=1)\n",
    "    data['alt_pred'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].mean(axis=1)\n",
    "\n",
    "    data['pref_allele'] = data[['fdrp_bh_ref',\n",
    "           'fdrp_bh_alt']].apply(lambda x: 'ref' if x.fdrp_bh_ref < x.fdrp_bh_alt else \"alt\", axis=1)\n",
    "\n",
    "    data['score'] = np.log2(data['min_fdr'] + eps) * (data['pref_allele'].apply(lambda x: -1 if x=='ref' else 1))\n",
    "\n",
    "\n",
    "    var_ref = data.loc[:, data.columns.str.startswith(\"pred_ref\")].var(axis=1).mean()\n",
    "    var_alt = data.loc[:, data.columns.str.startswith(\"pred_alt\")].var(axis=1).mean()\n",
    "    var = var_ref + var_alt\n",
    "    std = np.sqrt(var / N)\n",
    "    #se = std * 1.64485\n",
    "    \n",
    "    data['pred_score'] = data['ref_pred'] - data['alt_pred']\n",
    "    \n",
    "    data['pred_pval'] = norm.cdf(-abs(data['pred_score']) / std)\n",
    "\n",
    "\n",
    "    print(\"All data, Pearson R2\", pearsonr(data['score'], data['ref_pred'] - data['alt_pred'])[0])\n",
    "\n",
    "    significant = abs(data['ref_pred'] - data['alt_pred']) > se \n",
    "    data['pred_pref'] = (data['ref_pred'] - data['alt_pred']).apply(lambda x: \"ref\" if x > 0 else \"alt\")\n",
    "    \n",
    "    mask = data['min_fdr'] < 0.05\n",
    "    \n",
    "    print(\"FDR < 0.05, data, Pearson R2\", pearsonr(data[mask]['score'], data[mask]['ref_pred'] - data[mask]['alt_pred'])[0])\n",
    "    \n",
    "    \n",
    "    data_sig = data[mask].copy()\n",
    "    crosstab = pd.crosstab(index=data['pref_allele'], columns=data['pred_pref_allele'])\n",
    "    crosstab = crosstab.loc[['ref', 'alt'], ['ref', 'alt']]\n",
    "    print(\"Data without thresholding model predictions\")\n",
    "    print(crosstab)\n",
    "    \n",
    "    \n",
    "    data.loc[~significant, 'pred_pref'] = \"no\"\n",
    "    data.loc[data['min_fdr'] > 0.05, 'pref_allele'] = \"no\"\n",
    "    crosstab = pd.crosstab(index=data['pref_allele'], columns=data['pred_pref'])\n",
    "    crosstab = crosstab.loc[['ref', 'alt', 'no'], ['ref', 'alt', 'no']]\n",
    "    print(crosstab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "361fe312",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f41b5c07",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "1d913d14",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All data PearsonRResult(statistic=0.07868945961330985, pvalue=0.0)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3579166/3563733882.py:28: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  data['pred_pref'][~significant] = \"no\"\n",
      "/tmp/ipykernel_3579166/3563733882.py:29: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  data['pref_allele'][data['min_fdr'] > 0.05] = \"no\"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Signigicant data PearsonRResult(statistic=0.25821531396716446, pvalue=9.680540476622922e-146)\n"
     ]
    },
    {
     "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>pred_pref</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>no</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pref_allele</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ref</th>\n",
       "      <td>1944</td>\n",
       "      <td>1148</td>\n",
       "      <td>1726</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alt</th>\n",
       "      <td>1150</td>\n",
       "      <td>1936</td>\n",
       "      <td>1676</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>no</th>\n",
       "      <td>83457</td>\n",
       "      <td>83640</td>\n",
       "      <td>109619</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "pred_pref      ref    alt      no\n",
       "pref_allele                      \n",
       "ref           1944   1148    1726\n",
       "alt           1150   1936    1676\n",
       "no           83457  83640  109619"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eps =1e-10\n",
    "tp = \"single\"\n",
    "data = pd.read_csv(\"../human_legnet/HepG2_asb.tsv\", sep=\"\\t\")\n",
    "data['ref_pred'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].mean(axis=1)\n",
    "data['alt_pred'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].mean(axis=1)\n",
    "data['ref_pred_std'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].std(axis=1)\n",
    "data['alt_pred_std'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].std(axis=1)\n",
    "data['pref_allele'] = data[['es_mean_ref',\n",
    "       'es_mean_alt']].apply(lambda x: 'ref' if x.es_mean_ref > x.es_mean_alt else \"alt\", axis=1)\n",
    "\n",
    "data['score'] = np.log2(data['min_fdr'] + eps) * (data['pref_allele'].apply(lambda x: -1 if x=='ref' else 1))\n",
    "\n",
    "data['pred_std'] =  np.sqrt((data['ref_pred_std'] ** 2 + data['alt_pred_std'] ** 2) / 2)\n",
    "\n",
    "if tp == \"shared\":\n",
    "    std =  (data['ref_pred'] - data['alt_pred']).std()\n",
    "elif tp == \"single\":\n",
    "    std = np.sqrt(2) * data['pred_std'] / np.sqrt(data.columns.str.startswith(\"pred_alt\").sum())\n",
    "else:\n",
    "    raise NotImplementedError()\n",
    "\n",
    "data['pred_std'] =  np.sqrt((data['ref_pred_std'] ** 2 + data['alt_pred_std'] ** 2) / 2)\n",
    "\n",
    "print(\"All data\", pearsonr(data['score'], data['ref_pred'] - data['alt_pred']))\n",
    "\n",
    "significant = abs(data['ref_pred'] - data['alt_pred']) > std \n",
    "data['pred_pref'] = (data['ref_pred'] - data['alt_pred']).apply(lambda x: \"ref\" if x > 0 else \"alt\")\n",
    "data['pred_pref'][~significant] = \"no\"\n",
    "data['pref_allele'][data['min_fdr'] > 0.05] = \"no\"\n",
    "\n",
    "mask = data['min_fdr'] < 0.05\n",
    "\n",
    "print(\"Signigicant data\", pearsonr(data[mask]['score'], data[mask]['ref_pred'] - data[mask]['alt_pred']))\n",
    "crosstab = pd.crosstab(index=data['pref_allele'], columns=data['pred_pref'])\n",
    "crosstab = crosstab.loc[['ref', 'alt', 'no'], ['ref', 'alt', 'no']]\n",
    "crosstab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "2aa97ecb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All data PearsonRResult(statistic=0.07020284019780694, pvalue=0.0)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_3579166/2987647316.py:28: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  data['pred_pref'][~significant] = \"no\"\n",
      "/tmp/ipykernel_3579166/2987647316.py:29: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  data['pref_allele'][data['min_fdr'] > 0.05] = \"no\"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Signigicant data PearsonRResult(statistic=0.16583537426025294, pvalue=1.3864612161920364e-168)\n"
     ]
    },
    {
     "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>pred_pref</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>no</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>pref_allele</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ref</th>\n",
       "      <td>5325</td>\n",
       "      <td>3796</td>\n",
       "      <td>3278</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alt</th>\n",
       "      <td>4816</td>\n",
       "      <td>6172</td>\n",
       "      <td>4081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>no</th>\n",
       "      <td>104998</td>\n",
       "      <td>104093</td>\n",
       "      <td>91380</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "pred_pref       ref     alt     no\n",
       "pref_allele                       \n",
       "ref            5325    3796   3278\n",
       "alt            4816    6172   4081\n",
       "no           104998  104093  91380"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eps =1e-10\n",
    "tp = \"single\"\n",
    "data = pd.read_csv(\"../human_legnet/K562_asb.tsv\", sep=\"\\t\")\n",
    "data['ref_pred'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].mean(axis=1)\n",
    "data['alt_pred'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].mean(axis=1)\n",
    "data['ref_pred_std'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].std(axis=1)\n",
    "data['alt_pred_std'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].std(axis=1)\n",
    "data['pref_allele'] = data[['es_mean_ref',\n",
    "       'es_mean_alt']].apply(lambda x: 'ref' if x.es_mean_ref > x.es_mean_alt else \"alt\", axis=1)\n",
    "\n",
    "data['score'] = np.log2(data['min_fdr'] + eps) * (data['pref_allele'].apply(lambda x: -1 if x=='ref' else 1))\n",
    "\n",
    "data['pred_std'] =  np.sqrt((data['ref_pred_std'] ** 2 + data['alt_pred_std'] ** 2) / 2)\n",
    "\n",
    "if tp == \"shared\":\n",
    "    std =  (data['ref_pred'] - data['alt_pred']).std()\n",
    "elif tp == \"single\":\n",
    "    std = np.sqrt(2) * data['pred_std'] / np.sqrt(data.columns.str.startswith(\"pred_alt\").sum())\n",
    "else:\n",
    "    raise NotImplementedError()\n",
    "\n",
    "data['pred_std'] =  np.sqrt((data['ref_pred_std'] ** 2 + data['alt_pred_std'] ** 2) / 2)\n",
    "\n",
    "print(\"All data\", pearsonr(data['score'], data['ref_pred'] - data['alt_pred']))\n",
    "\n",
    "significant = abs(data['ref_pred'] - data['alt_pred']) > std \n",
    "data['pred_pref'] = (data['ref_pred'] - data['alt_pred']).apply(lambda x: \"ref\" if x > 0 else \"alt\")\n",
    "data['pred_pref'][~significant] = \"no\"\n",
    "data['pref_allele'][data['min_fdr'] > 0.05] = \"no\"\n",
    "\n",
    "mask = data['min_fdr'] < 0.05\n",
    "\n",
    "print(\"Signigicant data\", pearsonr(data[mask]['score'], data[mask]['ref_pred'] - data[mask]['alt_pred']))\n",
    "crosstab = pd.crosstab(index=data['pref_allele'], columns=data['pred_pref'])\n",
    "crosstab = crosstab.loc[['ref', 'alt', 'no'], ['ref', 'alt', 'no']]\n",
    "crosstab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "ae6a8f9c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "bf96e374",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['chr', 'pos', 'ID', 'ref', 'alt', 'repeat_type', 'n_peak_calls',\n",
       "       'n_peak_callers', 'mean_BAD', 'mean_SNP_per_segment', 'n_aggregated',\n",
       "       'total_cover', 'es_mean_ref', 'es_mean_alt', 'fdrp_bh_ref',\n",
       "       'fdrp_bh_alt', 'min_fdr', 'score', 'cls',\n",
       "       'pred_ref_test10_val6_0_forw'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.columns[0:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "b18e49ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "y_true = data['pref_allele'] != \"no\"\n",
    "y_pred = data['pred_pref'] != \"no\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "738007f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "       False       0.93      0.93      0.93     31259\n",
      "        True       0.14      0.14      0.14      2429\n",
      "\n",
      "    accuracy                           0.88     33688\n",
      "   macro avg       0.54      0.54      0.54     33688\n",
      "weighted avg       0.88      0.88      0.88     33688\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(classification_report(y_pred=y_pred, y_true=y_true))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "c606e987",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "772M\t../human_legnet/HepG2_asb.tsv\r\n",
      "47M\t../human_legnet/HepG2_cov_atac.tsv\r\n",
      "42M\t../human_legnet/HepG2_cov_dnase.tsv\r\n",
      "4,2M\t../human_legnet/HepG2_cov.tsv\r\n",
      "894M\t../human_legnet/K562_asb.tsv\r\n",
      "26M\t../human_legnet/K562_cov_atac.tsv\r\n",
      "96M\t../human_legnet/K562_cov_dnase.tsv\r\n",
      "234M\t../human_legnet/K562_cov.tsv\r\n",
      "40M\t../human_legnet/pred.tsv\r\n"
     ]
    }
   ],
   "source": [
    "!du -h -l ../human_legnet/*.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d84665b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tp = \"shared\"\n",
    "data = pd.read_csv(\"../human_legnet/HepG2_cov_dnase.tsv\", sep=\"\\t\")\n",
    "data['ref_pred'] = data.loc[:, data.columns.str.startswith(\"pred_ref\")].mean(axis=1)\n",
    "data['alt_pred'] = data.loc[:, data.columns.str.startswith(\"pred_alt\")].mean(axis=1)\n",
    "N =  data.columns.str.startswith(\"pred_ref\").sum()\n",
    "\n",
    "#data['pred_std'] =  np.sqrt((data['ref_pred_std'] ** 2 + data['alt_pred_std'] ** 2) / 2)\n",
    "\n",
    "if tp == \"shared\":\n",
    "    var_ref = data.loc[:, data.columns.str.startswith(\"pred_ref\")].var(axis=1).mean()\n",
    "    var_alt = data.loc[:, data.columns.str.startswith(\"pred_alt\")].var(axis=1).mean()\n",
    "    var = var_ref + var_alt\n",
    "    std = np.sqrt(var / N)\n",
    "    se = std * 1.64485\n",
    "elif tp == \"single\":\n",
    "    var_ref = data.loc[:, data.columns.str.startswith(\"pred_ref\")].var(axis=1)\n",
    "    var_alt = data.loc[:, data.columns.str.startswith(\"pred_alt\")].var(axis=1)\n",
    "    var = var_ref + var_alt\n",
    "    std = np.sqrt(var / N)\n",
    "    se = std * 1.64485\n",
    "else:\n",
    "    raise NotImplementedError()\n",
    "\n",
    "\n",
    "print(\"All data\", pearsonr(data['score'], data['ref_pred'] - data['alt_pred']))\n",
    "\n",
    "significant = abs(data['ref_pred'] - data['alt_pred']) > se \n",
    "data['pred_pref'] = (data['ref_pred'] - data['alt_pred']).apply(lambda x: \"ref\" if x > 0 else \"alt\")\n",
    "data['pred_pref'][~significant] = \"no\"\n",
    "data['pref_allele'][data['fdr_comb_pval'] > 0.05] = \"no\"\n",
    "\n",
    "mask = data['fdr_comb_pval'] < 0.05\n",
    "\n",
    "print(\"Signigicant data\", pearsonr(data[mask]['score'], data[mask]['ref_pred'] - data[mask]['alt_pred']))\n",
    "crosstab = pd.crosstab(index=data['pref_allele'], columns=data['pred_pref'])\n",
    "crosstab = crosstab.loc[['ref', 'alt', 'no'], ['ref', 'alt', 'no']]\n",
    "crosstab"
   ]
  }
 ],
 "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}