6_saving_and_basic_qc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os\n",
    "import glob\n",
    "import importlib\n",
    "import log_parser\n",
    "import scipy.stats as st\n",
    "import importlib\n",
    "\n",
    "from basic_tools import *\n",
    "\n",
    "debug=False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Loading h2 results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2_neale_result=pd.read_csv(h2_path,sep='\\t')\n",
    "h2_neale_result=h2_neale_result.set_index('phenotype')\n",
    "h2_neale_sorted=h2_neale_result.sort_values('h2_observed',ascending=False)\n",
    "h2_neale_sorted.columns = h2_neale_sorted.columns.map(lambda x: 'neale.'+str(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2_neale_v2_result=pd.read_csv(h2_v2_path,sep='\\t')\n",
    "h2_neale_v2_result=h2_neale_v2_result.set_index('phenotype')\n",
    "h2_neale_v2_sorted=h2_neale_v2_result.sort_values('h2_observed',ascending=False)\n",
    "h2_neale_v2_sorted.columns = h2_neale_v2_sorted.columns.map(lambda x: 'neale2.'+str(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## from my ldsc result of univariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>uni.filename</th>\n",
       "      <th>uni.Total Observed scale h2</th>\n",
       "      <th>uni.Total Observed scale h2_se</th>\n",
       "      <th>uni.Total Liability scale h2</th>\n",
       "      <th>uni.Total Liability scale h2_se</th>\n",
       "      <th>uni.Lambda GC</th>\n",
       "      <th>uni.Mean Chi^2</th>\n",
       "      <th>uni.Intercept</th>\n",
       "      <th>uni.Intercept_se</th>\n",
       "      <th>uni.Ratio</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>phenotype</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>50_irnt</th>\n",
       "      <td>data/out_final/uni.50_irnt.ldsc.log</td>\n",
       "      <td>0.3772</td>\n",
       "      <td>0.0197</td>\n",
       "      <td>0.3772</td>\n",
       "      <td>0.0197</td>\n",
       "      <td>2.5074</td>\n",
       "      <td>4.7497</td>\n",
       "      <td>1.6108</td>\n",
       "      <td>0.0373</td>\n",
       "      <td>0.1629 (0.0099)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20153_irnt</th>\n",
       "      <td>data/out_final/uni.20153_irnt.ldsc.log</td>\n",
       "      <td>0.3293</td>\n",
       "      <td>0.0209</td>\n",
       "      <td>0.3293</td>\n",
       "      <td>0.0209</td>\n",
       "      <td>1.5696</td>\n",
       "      <td>2.1139</td>\n",
       "      <td>1.2127</td>\n",
       "      <td>0.0188</td>\n",
       "      <td>0.1909 (0.0169)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5097_irnt</th>\n",
       "      <td>data/out_final/uni.5097_irnt.ldsc.log</td>\n",
       "      <td>0.2878</td>\n",
       "      <td>0.0233</td>\n",
       "      <td>0.2878</td>\n",
       "      <td>0.0233</td>\n",
       "      <td>1.3169</td>\n",
       "      <td>1.5435</td>\n",
       "      <td>1.1186</td>\n",
       "      <td>0.0134</td>\n",
       "      <td>0.2183 (0.0246)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5098_irnt</th>\n",
       "      <td>data/out_final/uni.5098_irnt.ldsc.log</td>\n",
       "      <td>0.2825</td>\n",
       "      <td>0.0223</td>\n",
       "      <td>0.2825</td>\n",
       "      <td>0.0223</td>\n",
       "      <td>1.3203</td>\n",
       "      <td>1.5455</td>\n",
       "      <td>1.1236</td>\n",
       "      <td>0.0135</td>\n",
       "      <td>0.2266 (0.0247)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5133_irnt</th>\n",
       "      <td>data/out_final/uni.5133_irnt.ldsc.log</td>\n",
       "      <td>0.2771</td>\n",
       "      <td>0.0221</td>\n",
       "      <td>0.2771</td>\n",
       "      <td>0.0221</td>\n",
       "      <td>1.3068</td>\n",
       "      <td>1.5356</td>\n",
       "      <td>1.1237</td>\n",
       "      <td>0.0135</td>\n",
       "      <td>0.231 (0.0252)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                      uni.filename  \\\n",
       "phenotype                                            \n",
       "50_irnt        data/out_final/uni.50_irnt.ldsc.log   \n",
       "20153_irnt  data/out_final/uni.20153_irnt.ldsc.log   \n",
       "5097_irnt    data/out_final/uni.5097_irnt.ldsc.log   \n",
       "5098_irnt    data/out_final/uni.5098_irnt.ldsc.log   \n",
       "5133_irnt    data/out_final/uni.5133_irnt.ldsc.log   \n",
       "\n",
       "            uni.Total Observed scale h2  uni.Total Observed scale h2_se  \\\n",
       "phenotype                                                                 \n",
       "50_irnt                          0.3772                          0.0197   \n",
       "20153_irnt                       0.3293                          0.0209   \n",
       "5097_irnt                        0.2878                          0.0233   \n",
       "5098_irnt                        0.2825                          0.0223   \n",
       "5133_irnt                        0.2771                          0.0221   \n",
       "\n",
       "            uni.Total Liability scale h2  uni.Total Liability scale h2_se  \\\n",
       "phenotype                                                                   \n",
       "50_irnt                           0.3772                           0.0197   \n",
       "20153_irnt                        0.3293                           0.0209   \n",
       "5097_irnt                         0.2878                           0.0233   \n",
       "5098_irnt                         0.2825                           0.0223   \n",
       "5133_irnt                         0.2771                           0.0221   \n",
       "\n",
       "            uni.Lambda GC  uni.Mean Chi^2  uni.Intercept  uni.Intercept_se  \\\n",
       "phenotype                                                                    \n",
       "50_irnt            2.5074          4.7497         1.6108            0.0373   \n",
       "20153_irnt         1.5696          2.1139         1.2127            0.0188   \n",
       "5097_irnt          1.3169          1.5435         1.1186            0.0134   \n",
       "5098_irnt          1.3203          1.5455         1.1236            0.0135   \n",
       "5133_irnt          1.3068          1.5356         1.1237            0.0135   \n",
       "\n",
       "                  uni.Ratio  \n",
       "phenotype                    \n",
       "50_irnt     0.1629 (0.0099)  \n",
       "20153_irnt  0.1909 (0.0169)  \n",
       "5097_irnt   0.2183 (0.0246)  \n",
       "5098_irnt   0.2266 (0.0247)  \n",
       "5133_irnt    0.231 (0.0252)  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h2_myuni_sorted=read_ldsc('uni',False)\n",
    "h2_myuni_sorted.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## from my ldsc result of various BP scales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bp300 was loaded\n",
      "bp128 was loaded\n",
      "bp64 was loaded\n",
      "bp32 was loaded\n",
      "bp16 was loaded\n",
      "bp8 was loaded\n"
     ]
    }
   ],
   "source": [
    "suffix_bp_list_power=\"bp300 bp128 bp64 bp32 bp16 bp8\".split()\n",
    "suffix_bp_list_full=\"bp300 bp128 bp64 bp32 bp16 bp8\".split()\n",
    "\n",
    "h2_bp_sorted_list=[]\n",
    "for suffix in suffix_bp_list_full:\n",
    "    h2_bp_sorted_list.append(read_ldsc(suffix,verbose=True))\n",
    "    print(suffix,\"was loaded\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## merging data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2_total=h2_bp_sorted_list[0].copy()\n",
    "for h2_bp_sorted in h2_bp_sorted_list[1:]:\n",
    "    h2_total=h2_total.join(h2_bp_sorted)\n",
    "\n",
    "h2_total=h2_total.join(h2_myuni_sorted)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2_total['uni.Total h2 z']=h2_total['uni.Total Observed scale h2']/h2_total['uni.Total Observed scale h2_se']\n",
    "h2_total['uni.Total h2 p']=st.norm.sf(h2_total['uni.Total Observed scale h2']/h2_total['uni.Total Observed scale h2_se'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2_total.loc[phenotypes_par_filtered.index].to_csv(h2_total_par_filtered_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'h2_neale_v2_sorted' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-1ec350e13156>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mh2_total\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mh2_total\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh2_neale_sorted\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mh2_total\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mh2_total\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjoin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh2_neale_v2_sorted\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'h2_neale_v2_sorted' is not defined"
     ]
    }
   ],
   "source": [
    "h2_total=h2_total.join(h2_neale_sorted)\n",
    "h2_total=h2_total.join(h2_neale_v2_sorted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Printing phenotypes that the three sources have in common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "186\n"
     ]
    }
   ],
   "source": [
    "#h2_total['neale.h2_observed'], h2_total['atlas.h2'])al['my.Total Observed scale h2'])\n",
    "conditions=[]\n",
    "conditions.append(np.isnan(h2_total['neale.h2_observed']))\n",
    "conditions.append(np.isnan(h2_total['neale2.h2_observed']))\n",
    "conditions.append(np.isnan(h2_total['uni.Total Observed scale h2']))\n",
    "\n",
    "for suffix in suffix_bp_list_full:\n",
    "    conditions.append(np.isnan(h2_total['{}.Total Observed scale h2'.format(suffix)]))    \n",
    "    \n",
    "common=np.logical_not(np.logical_or.reduce(conditions))\n",
    "print(np.sum(common))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "h2_total=h2_total.loc[phenotypes_par_filtered.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Verification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parsing result of univariate LDSC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finished loading uni\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'h2': 0.1718,\n",
       " 'h2_se': 0.0058,\n",
       " 'intercept': 1.2555,\n",
       " 'category': 'uni',\n",
       " 'M': 6856128}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pheno_code='23115_irnt'\n",
    "regression_result_uni=parse_uni_regression_result(h2_total,pheno_code)\n",
    "regression_result_uni"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parsing result of partitioned LDSC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "suffix_bp_list_power=\"bp300 bp128 bp64 bp32 bp16 bp8\".split()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "regression_result_bp_list_dict=dict([make_regression_result_list(h2_total,phenotype,suffix_bp_list_power) for phenotype in phenotypes_par_filtered.index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bp300: 0.1881\n",
      "bp128: 0.1895\n",
      "bp64: 0.19039999999999999\n",
      "bp32: 0.18999999999999997\n",
      "bp16: 0.18877615\n",
      "bp8: 0.18934133916999998\n"
     ]
    }
   ],
   "source": [
    "for idx,suffix in enumerate(suffix_bp_list_power):\n",
    "    print(\"{}: {}\".format(suffix,regression_result_bp_list_dict['23115_irnt'][idx]['h2'].sum()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "497"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(regression_result_bp_list_dict.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'h2_total_filtered' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-26-8002473d83a0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mh2_total_filtered\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mh2_total\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'h2_total_filtered' is not defined"
     ]
    }
   ],
   "source": [
    "h2_total_filtered.shape,h2_total.shape"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}