Lecture-10-SNPHeritability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estimating SNP heritability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup the environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from warnings import simplefilter \n",
    "simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "import sys\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.stats as st\n",
    "import pylab as pl\n",
    "import pandas as pd\n",
    "import h5py\n",
    "sp.random.seed(0)\n",
    "\n",
    "from limix.vardec import VarianceDecomposition\n",
    "import limix.util as lmx_util"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Load the data\n",
    "Load the arabidopsis data, which are in an hdf5 file.\n",
    "\n",
    "To process your own data, use the limix command line binary (see [here](http://nbviewer.jupyter.org/github/limix/limix-tutorials/blob/master/preprocessing_QC/loading_files.ipynb) for an example)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('./..')\n",
    "import data as tutorial_data\n",
    "file_name = tutorial_data.get_file('arab107')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up the data object\n",
    "The HDF5 file holds both the genotype and phenotype data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = h5py.File(file_name, 'r')\n",
    "phenotype_names = ['5_FT10','6_FT16','7_FT22']\n",
    "pheno_group = f['phenotype']\n",
    "pheno_df = pd.DataFrame(pheno_group['matrix'][:], \n",
    "                        columns=np.char.decode(pheno_group['col_header']['phenotype_ID'][:]),\n",
    "                        index=pheno_group['row_header']['sample_ID'][:])[phenotype_names]"
   ]
  },
  {
   "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>5_FT10</th>\n",
       "      <th>6_FT16</th>\n",
       "      <th>7_FT22</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>5837</th>\n",
       "      <td>57.0</td>\n",
       "      <td>50.0</td>\n",
       "      <td>43.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6008</th>\n",
       "      <td>60.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>24.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6009</th>\n",
       "      <td>98.0</td>\n",
       "      <td>151.0</td>\n",
       "      <td>250.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6016</th>\n",
       "      <td>75.0</td>\n",
       "      <td>101.0</td>\n",
       "      <td>113.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6040</th>\n",
       "      <td>71.0</td>\n",
       "      <td>92.0</td>\n",
       "      <td>87.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      5_FT10  6_FT16  7_FT22\n",
       "5837    57.0    50.0    43.0\n",
       "6008    60.0    41.0    24.0\n",
       "6009    98.0   151.0   250.0\n",
       "6016    75.0   101.0   113.0\n",
       "6040    71.0    92.0    87.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pheno_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>5_FT10</th>\n",
       "      <th>6_FT16</th>\n",
       "      <th>7_FT22</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>194.000000</td>\n",
       "      <td>193.000000</td>\n",
       "      <td>193.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>63.969072</td>\n",
       "      <td>64.647237</td>\n",
       "      <td>74.719689</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>17.821556</td>\n",
       "      <td>40.496864</td>\n",
       "      <td>71.747296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>41.000000</td>\n",
       "      <td>26.000000</td>\n",
       "      <td>23.300000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>49.000000</td>\n",
       "      <td>41.000000</td>\n",
       "      <td>30.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>59.250000</td>\n",
       "      <td>50.250000</td>\n",
       "      <td>44.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>71.000000</td>\n",
       "      <td>70.000000</td>\n",
       "      <td>75.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>121.000000</td>\n",
       "      <td>252.000000</td>\n",
       "      <td>250.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           5_FT10      6_FT16      7_FT22\n",
       "count  194.000000  193.000000  193.000000\n",
       "mean    63.969072   64.647237   74.719689\n",
       "std     17.821556   40.496864   71.747296\n",
       "min     41.000000   26.000000   23.300000\n",
       "25%     49.000000   41.000000   30.000000\n",
       "50%     59.250000   50.250000   44.000000\n",
       "75%     71.000000   70.000000   75.000000\n",
       "max    121.000000  252.000000  250.000000"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pheno_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "geno_group = f['genotype']\n",
    "chromosomes = geno_group['col_header']['chrom'][::10]\n",
    "positions = geno_group['col_header']['pos'][::10]\n",
    "geno_df = pd.DataFrame(geno_group['matrix'][:,::10], columns=positions, \n",
    "                       index=geno_group['row_header']['sample_ID'][:],\n",
    "                      dtype='float64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>657</th>\n",
       "      <th>7601</th>\n",
       "      <th>13045</th>\n",
       "      <th>25365</th>\n",
       "      <th>31926</th>\n",
       "      <th>41427</th>\n",
       "      <th>48118</th>\n",
       "      <th>55684</th>\n",
       "      <th>62259</th>\n",
       "      <th>69311</th>\n",
       "      <th>...</th>\n",
       "      <th>26931730</th>\n",
       "      <th>26934779</th>\n",
       "      <th>26937502</th>\n",
       "      <th>26941289</th>\n",
       "      <th>26945443</th>\n",
       "      <th>26951228</th>\n",
       "      <th>26957649</th>\n",
       "      <th>26964341</th>\n",
       "      <th>26969880</th>\n",
       "      <th>26973598</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>6122</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6121</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6119</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6116</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6115</th>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 21456 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      657       7601      13045     25365     31926     41427     48118     \\\n",
       "6122       0.0       0.0       0.0       0.0       1.0       0.0       0.0   \n",
       "6121       1.0       0.0       0.0       0.0       0.0       1.0       1.0   \n",
       "6119       0.0       0.0       0.0       0.0       0.0       0.0       0.0   \n",
       "6116       1.0       0.0       0.0       0.0       0.0       1.0       1.0   \n",
       "6115       1.0       0.0       0.0       0.0       0.0       1.0       1.0   \n",
       "\n",
       "      55684     62259     69311     ...  26931730  26934779  26937502  \\\n",
       "6122       1.0       0.0       0.0  ...       0.0       0.0       1.0   \n",
       "6121       0.0       0.0       1.0  ...       1.0       0.0       0.0   \n",
       "6119       1.0       0.0       0.0  ...       0.0       1.0       0.0   \n",
       "6116       0.0       0.0       1.0  ...       0.0       0.0       0.0   \n",
       "6115       0.0       0.0       1.0  ...       0.0       0.0       1.0   \n",
       "\n",
       "      26941289  26945443  26951228  26957649  26964341  26969880  26973598  \n",
       "6122       0.0       1.0       0.0       0.0       0.0       1.0       0.0  \n",
       "6121       0.0       1.0       0.0       0.0       0.0       0.0       0.0  \n",
       "6119       0.0       0.0       0.0       1.0       0.0       0.0       0.0  \n",
       "6116       0.0       0.0       1.0       1.0       0.0       0.0       1.0  \n",
       "6115       0.0       1.0       0.0       0.0       0.0       1.0       0.0  \n",
       "\n",
       "[5 rows x 21456 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geno_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 192 entries, 5837 to 100000\n",
      "Data columns (total 3 columns):\n",
      " #   Column  Non-Null Count  Dtype  \n",
      "---  ------  --------------  -----  \n",
      " 0   5_FT10  192 non-null    float64\n",
      " 1   6_FT16  192 non-null    float64\n",
      " 2   7_FT22  192 non-null    float64\n",
      "dtypes: float64(3)\n",
      "memory usage: 6.0 KB\n"
     ]
    }
   ],
   "source": [
    "#flowering phenotypes\n",
    "phenotype_ID =  pheno_df.columns[1]\n",
    "\n",
    "filtered_pheno_df = pheno_df.dropna()\n",
    "filtered_pheno_df.info()\n",
    "sample_idx = geno_df.index.intersection(filtered_pheno_df.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loaded 192 samples, 3 phenotypes, 21456 snps\n"
     ]
    }
   ],
   "source": [
    "phenotype_std = lmx_util.preprocess.rankStandardizeNormal(filtered_pheno_df.loc[sample_idx].values)\n",
    "snps = geno_df.loc[sample_idx]\n",
    "N = snps.shape[0]\n",
    "S = snps.shape[1]\n",
    "P = pheno_df.shape[1]\n",
    "print(\"loaded %d samples, %d phenotypes, %s snps\" % (N, P, S))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use the SNPs to estimate sample similarity/relatedness (also known as kinship)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [00:00<00:00, 2117.57it/s]\n"
     ]
    }
   ],
   "source": [
    "from limix.stats import linear_kinship, gower_norm\n",
    "sample_relatedness_unnormalized = linear_kinship(snps.values)\n",
    "sample_relatedness = gower_norm(sample_relatedness_unnormalized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.08575548, -0.00220991,  0.01914943, -0.0065966 , -0.02605002],\n",
       "       [-0.00220991,  0.98639869, -0.01216375, -0.01691818, -0.01916793],\n",
       "       [ 0.01914943, -0.01216375,  1.04626177, -0.0221134 ,  0.01371893],\n",
       "       [-0.0065966 , -0.01691818, -0.0221134 ,  1.31148388, -0.01127501],\n",
       "       [-0.02605002, -0.01916793,  0.01371893, -0.01127501,  0.93731567]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_relatedness[0:5, 0:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimating heritability\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimating the pseudo/narrow-sense heritability of: 5_FT10\n",
      "genetic: 0.90, noise: 0.10\n",
      "Estimating the pseudo/narrow-sense heritability of: 6_FT16\n",
      "genetic: 0.97, noise: 0.03\n",
      "Estimating the pseudo/narrow-sense heritability of: 7_FT22\n",
      "genetic: 1.00, noise: 0.00\n"
     ]
    }
   ],
   "source": [
    "for idx in range(len(phenotype_names)):\n",
    "    print(\"Estimating the pseudo/narrow-sense heritability of: \" + phenotype_names[idx])\n",
    "    y = phenotype_std[:,idx]\n",
    "    vc = VarianceDecomposition(y)\n",
    "    vc.addFixedEffect()\n",
    "    vc.addRandomEffect(K=sample_relatedness)\n",
    "    vc.addRandomEffect(is_noise=True)\n",
    "    vc.optimize()\n",
    "\n",
    "    _var = vc.getVarianceComps()\n",
    "    _var /= _var.sum()\n",
    "    _var = _var.ravel()\n",
    "    print('genetic: %.2f, noise: %.2f' % (_var[0],_var[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although flowering time is believed to be highly heritable, these values (90 - 100%) are likely to be overestimates. Note that measures of SNP-heritability also act as a measure of confounding due to population structure."
   ]
  }
 ],
 "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}