3.2_making_top_snp_table_META.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e8bbd76e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "from collections import defaultdict\n",
    "from tqdm import tqdm\n",
    "from pprint import pprint\n",
    "\n",
    "from Bio import Entrez\n",
    "Entrez.email = \"sample@bioinf.me\"\n",
    "\n",
    "from functools import lru_cache\n",
    "\n",
    "pd.options.mode.chained_assignment = None  # default='warn'\n",
    "tqdm.pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b9636f",
   "metadata": {},
   "source": [
    "## For meta data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9ff98669",
   "metadata": {},
   "outputs": [],
   "source": [
    "@lru_cache(maxsize=4096)\n",
    "def return_gene_by_id(uid):\n",
    "    handle = Entrez.esummary(db=\"gene\", id=uid)\n",
    "    uid_record = Entrez.read(handle)\n",
    "    handle.close()\n",
    "    uid_summary = uid_record[\"DocumentSummarySet\"][\"DocumentSummary\"][0]\n",
    "    return uid_summary[\"Name\"]\n",
    "\n",
    "\n",
    "@lru_cache(maxsize=4096)\n",
    "def return_gene_by_rsid(snp_id):\n",
    "    answer = []\n",
    "    record = Entrez.read(\n",
    "        Entrez.elink(dbfrom=\"snp\", id=snp_id.replace(\"rs\", \"\"), db=\"gene\")\n",
    "    )\n",
    "    results = record[0][\"LinkSetDb\"]\n",
    "    if len(results) < 1:\n",
    "        return None\n",
    "    results = results[0][\"Link\"]\n",
    "    for result in results:\n",
    "        uid = result[\"Id\"]\n",
    "        gene_name = return_gene_by_id(uid)\n",
    "        answer.append(gene_name)\n",
    "    return \",\".join(answer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "43e9b49f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/achangalidi/miniconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3444: DtypeWarning: Columns (1) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "/home/achangalidi/miniconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3444: DtypeWarning: Columns (5) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "/home/achangalidi/miniconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3444: DtypeWarning: Columns (0) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "100%|████████████████████████████████████████████████████████████████████████████████████████████████| 7/7 [00:05<00:00,  1.33it/s]\n",
      "100%|████████████████████████████████████████████████████████████████████████████████████████████| 103/103 [00:55<00:00,  1.86it/s]\n",
      "100%|████████████████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:01<00:00,  1.89it/s]\n"
     ]
    }
   ],
   "source": [
    "cutoff = 0.05\n",
    "top_features = [\"I9_HYPTENSPREG\", \"GEST_DIABETES\", \"O15_PRETERM\"]\n",
    "\n",
    "datas = []\n",
    "all_datas = defaultdict()\n",
    "ns = defaultdict()\n",
    "\n",
    "for feat in top_features:\n",
    "\n",
    "    cur_dir = \"analysis_n_\" + feat\n",
    "\n",
    "    cur_files = [f for f in os.listdir(f\"{cur_dir}/data\")]\n",
    "    metal_data_name = [f for f in cur_files if f[-4:] == \".TBL\" and \"extended_\" in f][0]\n",
    "    mtag_data_name = \"_mtag_meta.txt\"\n",
    "    ukb_data_name = [f for f in cur_files if \"gwas.imputed_v3\" in f][0]\n",
    "    finn_data_name = [f for f in cur_files if \"hg19lifted\" in f][0]\n",
    "\n",
    "    metal_data = pd.read_csv(f\"{cur_dir}/data/{metal_data_name}\", sep=\"\\t\").add_prefix(\n",
    "        \"metal_\"\n",
    "    )\n",
    "    mtag_data = (\n",
    "        pd.read_csv(f\"{cur_dir}/data/{mtag_data_name}\", sep=\"\\t\")\n",
    "        .rename(\n",
    "            columns={\n",
    "                \"SNP\": \"rsid\",\n",
    "                \"CHR\": \"chr\",\n",
    "                \"BP\": \"pos\",\n",
    "                \"A1\": \"ref\",\n",
    "                \"A2\": \"alt\",\n",
    "                \"meta_freq\": \"maf\",\n",
    "                \"mtag_beta\": \"beta\",\n",
    "                \"mtag_se\": \"se\",\n",
    "                \"mtag_z\": \"z\",\n",
    "                \"mtag_pval\": \"pval\",\n",
    "            }\n",
    "        )\n",
    "        .add_prefix(\"mtag_\")\n",
    "    )\n",
    "    ukb_data = pd.read_csv(f\"{cur_dir}/data/{ukb_data_name}\", sep=\"\\t\").add_prefix(\n",
    "        \"ukb_\"\n",
    "    )\n",
    "    finn_data = pd.read_csv(f\"{cur_dir}/data/{finn_data_name}\", sep=\"\\t\").add_prefix(\n",
    "        \"finn_\"\n",
    "    )\n",
    "\n",
    "    metal_data[\"metal_significant\"] = (\n",
    "        metal_data.metal_pval < cutoff / metal_data.shape[0]\n",
    "    )\n",
    "    mtag_data[\"mtag_significant\"] = mtag_data.mtag_pval < cutoff / mtag_data.shape[0]\n",
    "    ns[feat] = [\n",
    "        metal_data.shape[0],\n",
    "        mtag_data.shape[0],\n",
    "        ukb_data.shape[0],\n",
    "        finn_data.shape[0],\n",
    "    ]\n",
    "\n",
    "    data = pd.merge(\n",
    "        metal_data, mtag_data, how=\"outer\", left_on=\"metal_rsid\", right_on=\"mtag_rsid\"\n",
    "    )\n",
    "    all_datas[feat] = data\n",
    "\n",
    "    data = data[\n",
    "        data.mtag_significant\n",
    "        | data.metal_significant\n",
    "        | (data.mtag_significant.isna() & data.metal_significant)\n",
    "    ]\n",
    "    data = pd.merge(\n",
    "        data, ukb_data, how=\"left\", left_on=\"metal_rsid\", right_on=\"ukb_rsid\"\n",
    "    )\n",
    "    data = pd.merge(\n",
    "        data, finn_data, how=\"left\", left_on=\"metal_rsid\", right_on=\"finn_rsid\"\n",
    "    )\n",
    "    data = data.rename(\n",
    "        columns={\n",
    "            \"metal_rsid\": \"rsid\",\n",
    "            \"metal_chr\": \"chr\",\n",
    "            \"metal_pos\": \"pos\",\n",
    "            \"metal_ref\": \"ref\",\n",
    "            \"metal_alt\": \"alt\",\n",
    "        }\n",
    "    )\n",
    "    data[\"trait\"] = feat\n",
    "    data[\"gene\"] = data.rsid.progress_apply(return_gene_by_rsid)\n",
    "    data = data[\n",
    "        [\n",
    "            \"rsid\",\n",
    "            \"chr\",\n",
    "            \"gene\",\n",
    "            \"pos\",\n",
    "            \"ref\",\n",
    "            \"alt\",\n",
    "            \"metal_pval\",\n",
    "            \"mtag_pval\",\n",
    "            \"ukb_pval\",\n",
    "            \"finn_pval\",\n",
    "            \"trait\",\n",
    "            \"metal_significant\",\n",
    "            \"mtag_significant\",\n",
    "        ]\n",
    "    ]\n",
    "    datas.append(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e9d61f1",
   "metadata": {},
   "source": [
    "## Prepare to sort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f3e67a26",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'1': 0,\n",
       " '2': 1,\n",
       " '3': 2,\n",
       " '4': 3,\n",
       " '5': 4,\n",
       " '6': 5,\n",
       " '7': 6,\n",
       " '8': 7,\n",
       " '9': 8,\n",
       " '10': 9,\n",
       " '11': 10,\n",
       " '12': 11,\n",
       " '13': 12,\n",
       " '14': 13,\n",
       " '15': 14,\n",
       " '16': 15,\n",
       " '17': 16,\n",
       " '18': 17,\n",
       " '19': 18,\n",
       " '20': 19,\n",
       " '21': 20,\n",
       " '22': 21,\n",
       " 'X': 22}"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import re\n",
    "\n",
    "def atof(text):\n",
    "    try:\n",
    "        retval = float(text)\n",
    "    except ValueError:\n",
    "        retval = text\n",
    "    return retval\n",
    "\n",
    "def natural_keys(text):\n",
    "    '''\n",
    "    alist.sort(key=natural_keys) sorts in human order\n",
    "    http://nedbatchelder.com/blog/200712/human_sorting.html\n",
    "    (See Toothy's implementation in the comments)\n",
    "    float regex comes from https://stackoverflow.com/a/12643073/190597\n",
    "    '''\n",
    "    return [ atof(c) for c in re.split(r'[+-]?([0-9]+(?:[.][0-9]*)?|[.][0-9]+)', text) ]\n",
    "kk = list(all_datas[feat].metal_chr.astype(str).unique())\n",
    "kk.sort(key=natural_keys)\n",
    "order = {i:j for i, j in zip(kk, range(len(kk)))}\n",
    "order"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8bbe6a3",
   "metadata": {},
   "source": [
    "## Sort and polish"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5df9cc4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "STARTING_GROUP = 0\n",
    "\n",
    "def make_group_col(df, col, group_col='group'):\n",
    "    global STARTING_GROUP\n",
    "    df[group_col] = None\n",
    "    df.groupby((df[col].shift() != df[col]).cumsum())\n",
    "    for k, v in df.groupby((df[col].shift() != df[col]).cumsum()):\n",
    "        df.loc[v.index,group_col] = STARTING_GROUP+k\n",
    "    STARTING_GROUP += k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "173f5eda",
   "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>rsid</th>\n",
       "      <th>chr</th>\n",
       "      <th>gene</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>metal_pval</th>\n",
       "      <th>mtag_pval</th>\n",
       "      <th>ukb_pval</th>\n",
       "      <th>finn_pval</th>\n",
       "      <th>trait</th>\n",
       "      <th>metal_significant</th>\n",
       "      <th>mtag_significant</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rs10882398</td>\n",
       "      <td>10</td>\n",
       "      <td>PLCE1</td>\n",
       "      <td>95892788</td>\n",
       "      <td>T</td>\n",
       "      <td>A</td>\n",
       "      <td>8.943000e-09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.005538</td>\n",
       "      <td>3.480000e-09</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>rs10882397</td>\n",
       "      <td>10</td>\n",
       "      <td>PLCE1</td>\n",
       "      <td>95892659</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>1.456000e-08</td>\n",
       "      <td>7.309473e-09</td>\n",
       "      <td>0.003695</td>\n",
       "      <td>1.950000e-08</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>rs16998073</td>\n",
       "      <td>4</td>\n",
       "      <td>None</td>\n",
       "      <td>81184341</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>9.128000e-09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.035201</td>\n",
       "      <td>1.910000e-11</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>rs932764</td>\n",
       "      <td>10</td>\n",
       "      <td>PLCE1</td>\n",
       "      <td>95895940</td>\n",
       "      <td>A</td>\n",
       "      <td>G</td>\n",
       "      <td>1.388000e-08</td>\n",
       "      <td>7.024045e-09</td>\n",
       "      <td>0.002932</td>\n",
       "      <td>2.980000e-08</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>rs12509595</td>\n",
       "      <td>4</td>\n",
       "      <td>None</td>\n",
       "      <td>81182554</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>1.575000e-08</td>\n",
       "      <td>7.070137e-09</td>\n",
       "      <td>0.043569</td>\n",
       "      <td>2.570000e-11</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>rs35954793</td>\n",
       "      <td>4</td>\n",
       "      <td>FGF5</td>\n",
       "      <td>81188513</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>6.145000e-09</td>\n",
       "      <td>2.680313e-09</td>\n",
       "      <td>0.029762</td>\n",
       "      <td>1.600000e-11</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>rs167479</td>\n",
       "      <td>19</td>\n",
       "      <td>RGL3</td>\n",
       "      <td>11526765</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>5.233000e-09</td>\n",
       "      <td>2.223811e-09</td>\n",
       "      <td>0.038050</td>\n",
       "      <td>5.130000e-12</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         rsid chr   gene       pos ref alt    metal_pval     mtag_pval  \\\n",
       "0  rs10882398  10  PLCE1  95892788   T   A  8.943000e-09           NaN   \n",
       "1  rs10882397  10  PLCE1  95892659   C   A  1.456000e-08  7.309473e-09   \n",
       "2  rs16998073   4   None  81184341   A   T  9.128000e-09           NaN   \n",
       "3    rs932764  10  PLCE1  95895940   A   G  1.388000e-08  7.024045e-09   \n",
       "4  rs12509595   4   None  81182554   T   C  1.575000e-08  7.070137e-09   \n",
       "5  rs35954793   4   FGF5  81188513   C   A  6.145000e-09  2.680313e-09   \n",
       "6    rs167479  19   RGL3  11526765   G   T  5.233000e-09  2.223811e-09   \n",
       "\n",
       "   ukb_pval     finn_pval           trait  metal_significant mtag_significant  \n",
       "0  0.005538  3.480000e-09  I9_HYPTENSPREG               True              NaN  \n",
       "1  0.003695  1.950000e-08  I9_HYPTENSPREG              False             True  \n",
       "2  0.035201  1.910000e-11  I9_HYPTENSPREG               True              NaN  \n",
       "3  0.002932  2.980000e-08  I9_HYPTENSPREG              False             True  \n",
       "4  0.043569  2.570000e-11  I9_HYPTENSPREG              False             True  \n",
       "5  0.029762  1.600000e-11  I9_HYPTENSPREG               True             True  \n",
       "6  0.038050  5.130000e-12  I9_HYPTENSPREG               True             True  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datas[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "16be085a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(datas)):\n",
    "    datas[i].pos = datas[i].pos.astype(int)\n",
    "    datas[i].chr = datas[i].chr.astype(str)\n",
    "    datas[i]['chr_order'] = datas[i].chr.apply(lambda x: order[x])\n",
    "    datas[i]['chr_gene'] = datas[i].chr + '_' + datas[i].gene.astype(str)\n",
    "    datas[i] = datas[i].sort_values(by=['chr_order', 'pos'])\n",
    "    make_group_col(datas[i], 'chr_gene', 'group')\n",
    "    datas[i] = datas[i].drop(['chr_order', 'chr_gene'], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "076f92d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9, 14)\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></th>\n",
       "      <th>rsid</th>\n",
       "      <th>chr</th>\n",
       "      <th>gene</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>metal_pval</th>\n",
       "      <th>mtag_pval</th>\n",
       "      <th>ukb_pval</th>\n",
       "      <th>finn_pval</th>\n",
       "      <th>trait</th>\n",
       "      <th>metal_significant</th>\n",
       "      <th>mtag_significant</th>\n",
       "      <th>group</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>rs16998073</td>\n",
       "      <td>4</td>\n",
       "      <td>None</td>\n",
       "      <td>81184341</td>\n",
       "      <td>A</td>\n",
       "      <td>T</td>\n",
       "      <td>9.128000e-09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.035201</td>\n",
       "      <td>1.910000e-11</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>rs35954793</td>\n",
       "      <td>4</td>\n",
       "      <td>FGF5</td>\n",
       "      <td>81188513</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>6.145000e-09</td>\n",
       "      <td>2.680313e-09</td>\n",
       "      <td>0.029762</td>\n",
       "      <td>1.600000e-11</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rs10882398</td>\n",
       "      <td>10</td>\n",
       "      <td>PLCE1</td>\n",
       "      <td>95892788</td>\n",
       "      <td>T</td>\n",
       "      <td>A</td>\n",
       "      <td>8.943000e-09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.005538</td>\n",
       "      <td>3.480000e-09</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>rs167479</td>\n",
       "      <td>19</td>\n",
       "      <td>RGL3</td>\n",
       "      <td>11526765</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>5.233000e-09</td>\n",
       "      <td>2.223811e-09</td>\n",
       "      <td>0.038050</td>\n",
       "      <td>5.130000e-12</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>rs36090025</td>\n",
       "      <td>10</td>\n",
       "      <td>TCF7L2</td>\n",
       "      <td>114774433</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>3.483000e-15</td>\n",
       "      <td>2.715405e-15</td>\n",
       "      <td>0.002077</td>\n",
       "      <td>1.910000e-18</td>\n",
       "      <td>GEST_DIABETES</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>101</th>\n",
       "      <td>rs7945617</td>\n",
       "      <td>11</td>\n",
       "      <td>None</td>\n",
       "      <td>92700287</td>\n",
       "      <td>T</td>\n",
       "      <td>C</td>\n",
       "      <td>2.798000e-31</td>\n",
       "      <td>2.037698e-31</td>\n",
       "      <td>0.036522</td>\n",
       "      <td>1.110000e-56</td>\n",
       "      <td>GEST_DIABETES</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>rs10830963</td>\n",
       "      <td>11</td>\n",
       "      <td>MTNR1B</td>\n",
       "      <td>92708710</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>4.530000e-41</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.136495</td>\n",
       "      <td>7.860000e-84</td>\n",
       "      <td>GEST_DIABETES</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>rs138109677</td>\n",
       "      <td>11</td>\n",
       "      <td>None</td>\n",
       "      <td>92725992</td>\n",
       "      <td>A</td>\n",
       "      <td>AATGTT</td>\n",
       "      <td>7.213000e-15</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.440159</td>\n",
       "      <td>7.060000e-30</td>\n",
       "      <td>GEST_DIABETES</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rs2963457</td>\n",
       "      <td>5</td>\n",
       "      <td>None</td>\n",
       "      <td>157907974</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>6.520000e-09</td>\n",
       "      <td>4.715197e-09</td>\n",
       "      <td>0.000204</td>\n",
       "      <td>1.940000e-06</td>\n",
       "      <td>O15_PRETERM</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            rsid chr    gene        pos ref     alt    metal_pval  \\\n",
       "2     rs16998073   4    None   81184341   A       T  9.128000e-09   \n",
       "5     rs35954793   4    FGF5   81188513   C       A  6.145000e-09   \n",
       "0     rs10882398  10   PLCE1   95892788   T       A  8.943000e-09   \n",
       "6       rs167479  19    RGL3   11526765   G       T  5.233000e-09   \n",
       "33    rs36090025  10  TCF7L2  114774433   A       C  3.483000e-15   \n",
       "101    rs7945617  11    None   92700287   T       C  2.798000e-31   \n",
       "2     rs10830963  11  MTNR1B   92708710   C       G  4.530000e-41   \n",
       "4    rs138109677  11    None   92725992   A  AATGTT  7.213000e-15   \n",
       "0      rs2963457   5    None  157907974   C       T  6.520000e-09   \n",
       "\n",
       "        mtag_pval  ukb_pval     finn_pval           trait  metal_significant  \\\n",
       "2             NaN  0.035201  1.910000e-11  I9_HYPTENSPREG               True   \n",
       "5    2.680313e-09  0.029762  1.600000e-11  I9_HYPTENSPREG               True   \n",
       "0             NaN  0.005538  3.480000e-09  I9_HYPTENSPREG               True   \n",
       "6    2.223811e-09  0.038050  5.130000e-12  I9_HYPTENSPREG               True   \n",
       "33   2.715405e-15  0.002077  1.910000e-18   GEST_DIABETES               True   \n",
       "101  2.037698e-31  0.036522  1.110000e-56   GEST_DIABETES               True   \n",
       "2             NaN  0.136495  7.860000e-84   GEST_DIABETES               True   \n",
       "4             NaN  0.440159  7.060000e-30   GEST_DIABETES               True   \n",
       "0    4.715197e-09  0.000204  1.940000e-06     O15_PRETERM               True   \n",
       "\n",
       "    mtag_significant group  \n",
       "2                NaN     1  \n",
       "5               True     2  \n",
       "0                NaN     3  \n",
       "6               True     4  \n",
       "33              True     5  \n",
       "101             True     6  \n",
       "2                NaN     7  \n",
       "4                NaN     8  \n",
       "0               True     9  "
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_datas = []\n",
    "PVAL_COL = 'metal_pval'\n",
    "for data in datas:\n",
    "    data\n",
    "    best = data.groupby('group').agg(\n",
    "        pval=pd.NamedAgg(column=PVAL_COL, aggfunc=\"min\"))\n",
    "    filtered_data = data[data[PVAL_COL].isin(best.pval)]\n",
    "    # может произойти что в одной хромосоме данный пвал макс, а в другой нет \n",
    "    # - важно взять нужную и только ее! \n",
    "    flags = []\n",
    "    for key, value in filtered_data.iterrows():\n",
    "        flags.append(best[best.pval==value[PVAL_COL]].index[0]==value.group)\n",
    "    best_datas.append(filtered_data[flags])\n",
    "final = pd.concat(best_datas)\n",
    "final.to_csv('./data/meta_top.tsv', sep='\\t', index=False)\n",
    "print(final.shape)\n",
    "final"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63c5c7fe",
   "metadata": {},
   "source": [
    "Next, you should remove unnecessary SNPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "8a44dbbf",
   "metadata": {
    "scrolled": true
   },
   "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>rsid</th>\n",
       "      <th>chr</th>\n",
       "      <th>gene</th>\n",
       "      <th>pos</th>\n",
       "      <th>ref</th>\n",
       "      <th>alt</th>\n",
       "      <th>metal_pval</th>\n",
       "      <th>mtag_pval</th>\n",
       "      <th>ukb_pval</th>\n",
       "      <th>finn_pval</th>\n",
       "      <th>trait</th>\n",
       "      <th>metal_significant</th>\n",
       "      <th>mtag_significant</th>\n",
       "      <th>effect</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>rs35954793</td>\n",
       "      <td>4</td>\n",
       "      <td>FGF5</td>\n",
       "      <td>81188513</td>\n",
       "      <td>C</td>\n",
       "      <td>A</td>\n",
       "      <td>6.145000e-09</td>\n",
       "      <td>2.680313e-09</td>\n",
       "      <td>0.029762</td>\n",
       "      <td>1.600000e-11</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>FGF5 : Intron Variant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>rs10882398</td>\n",
       "      <td>10</td>\n",
       "      <td>PLCE1</td>\n",
       "      <td>95892788</td>\n",
       "      <td>T</td>\n",
       "      <td>A</td>\n",
       "      <td>8.943000e-09</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.005538</td>\n",
       "      <td>3.480000e-09</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>PLCE1 : Intron Variant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>rs167479</td>\n",
       "      <td>19</td>\n",
       "      <td>RGL3</td>\n",
       "      <td>11526765</td>\n",
       "      <td>G</td>\n",
       "      <td>T</td>\n",
       "      <td>5.233000e-09</td>\n",
       "      <td>2.223811e-09</td>\n",
       "      <td>0.038050</td>\n",
       "      <td>5.130000e-12</td>\n",
       "      <td>I9_HYPTENSPREG</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>RGL3 : Missense Variant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>rs36090025</td>\n",
       "      <td>10</td>\n",
       "      <td>TCF7L2</td>\n",
       "      <td>114774433</td>\n",
       "      <td>A</td>\n",
       "      <td>C</td>\n",
       "      <td>3.483000e-15</td>\n",
       "      <td>2.715405e-15</td>\n",
       "      <td>0.002077</td>\n",
       "      <td>1.910000e-18</td>\n",
       "      <td>GEST_DIABETES</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>TCF7L2 : Intron Variant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>rs10830963</td>\n",
       "      <td>11</td>\n",
       "      <td>MTNR1B</td>\n",
       "      <td>92708710</td>\n",
       "      <td>C</td>\n",
       "      <td>G</td>\n",
       "      <td>4.530000e-41</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.136495</td>\n",
       "      <td>7.860000e-84</td>\n",
       "      <td>GEST_DIABETES</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "      <td>MTNR1B : Intron Variant</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>rs2963457</td>\n",
       "      <td>5</td>\n",
       "      <td>EBF1</td>\n",
       "      <td>157907974</td>\n",
       "      <td>C</td>\n",
       "      <td>T</td>\n",
       "      <td>6.520000e-09</td>\n",
       "      <td>4.715197e-09</td>\n",
       "      <td>0.000204</td>\n",
       "      <td>1.940000e-06</td>\n",
       "      <td>O15_PRETERM</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         rsid  chr    gene        pos ref alt    metal_pval     mtag_pval  \\\n",
       "0  rs35954793    4    FGF5   81188513   C   A  6.145000e-09  2.680313e-09   \n",
       "1  rs10882398   10   PLCE1   95892788   T   A  8.943000e-09           NaN   \n",
       "2    rs167479   19    RGL3   11526765   G   T  5.233000e-09  2.223811e-09   \n",
       "3  rs36090025   10  TCF7L2  114774433   A   C  3.483000e-15  2.715405e-15   \n",
       "4  rs10830963   11  MTNR1B   92708710   C   G  4.530000e-41           NaN   \n",
       "5   rs2963457    5    EBF1  157907974   C   T  6.520000e-09  4.715197e-09   \n",
       "\n",
       "   ukb_pval     finn_pval           trait  metal_significant mtag_significant  \\\n",
       "0  0.029762  1.600000e-11  I9_HYPTENSPREG               True             True   \n",
       "1  0.005538  3.480000e-09  I9_HYPTENSPREG               True              NaN   \n",
       "2  0.038050  5.130000e-12  I9_HYPTENSPREG               True             True   \n",
       "3  0.002077  1.910000e-18   GEST_DIABETES               True             True   \n",
       "4  0.136495  7.860000e-84   GEST_DIABETES               True              NaN   \n",
       "5  0.000204  1.940000e-06     O15_PRETERM               True             True   \n",
       "\n",
       "                    effect  \n",
       "0    FGF5 : Intron Variant  \n",
       "1   PLCE1 : Intron Variant  \n",
       "2  RGL3 : Missense Variant  \n",
       "3  TCF7L2 : Intron Variant  \n",
       "4  MTNR1B : Intron Variant  \n",
       "5                      NaN  "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_final = pd.read_csv('./data/meta_top_short.csv', sep='\\t')\n",
    "_final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "d8ab0667",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c('rs35954793', 'rs10882398', 'rs167479'),\n",
      "c('rs36090025', 'rs10830963'),\n",
      "c('rs2963457')\n",
      "\n",
      "c('FGF5', 'PLCE1', 'RGL3'),\n",
      "c('TCF7L2', 'MTNR1B'),\n",
      "c('EBF1')\n",
      "\n",
      "'_MET_I9_HYPTENSPREG',\n",
      "'_MET_GEST_DIABETES',\n",
      "'_MET_O15_PRETERM'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "traits = _final.trait.unique()\n",
    "s=''\n",
    "for t in traits:\n",
    "    s+=\"c('\"+\"', '\".join(list(_final[_final.trait==t].rsid))+\"'),\\n\"\n",
    "print(s[:-2], end='\\n\\n')\n",
    "s=''\n",
    "for t in traits:\n",
    "    s+=\"c('\"+\"', '\".join(list(_final[_final.trait==t].gene))+\"'),\\n\"\n",
    "print(s[:-2], end='\\n\\n')\n",
    "s=''\n",
    "# for t in traits:\n",
    "#     s+=f\"'{DIR}/{next(filter(lambda f: t in f and 'flc_fmaf1' in f, files))}',\\n\"\n",
    "# print(s[:-2], end='\\n\\n')\n",
    "# s=''\n",
    "for t in traits:\n",
    "    s+=f\"'_MET_{t}',\\n\"\n",
    "print(s[:-2], end='\\n\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c10eff6",
   "metadata": {},
   "source": [
    "## Sort in chr:pos ascending order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "7c188094",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in top_features:\n",
    "    d = all_datas[t].rename(columns={'metal_rsid': 'rsid', 'metal_chr': 'chr', 'metal_pos': 'pos', 'metal_ref': 'ref', 'metal_alt': 'alt', 'metal_pval':'pval'})\n",
    "    d.pos = d.pos.astype(int)\n",
    "    d.chr = d.chr.astype(str)\n",
    "    d['chr_order'] = d.chr.apply(lambda x: order[x])\n",
    "    kkkk = d.sort_values(by=['chr_order', 'pos'])[['rsid', 'chr', 'pos', 'ref', 'alt', 'pval']]\n",
    "    cur_f = f\"meta_special/{t}.tsv\"\n",
    "    kkkk.to_csv(cur_f, sep='\\t', index=False)"
   ]
  }
 ],
 "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}