03-relatedness-estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "72b68b58",
   "metadata": {},
   "source": [
    "# Institute for Behavioral Genetics International Statistical Genetics 2023 Workshop \n",
    "\n",
    "## Genetic relatedness exploration with Hail\n",
    "\n",
    "In this practical, we will learn how to:\n",
    "\n",
    "1) Simulate random mating in Hail\n",
    "\n",
    "2) Run various relatedness estimation methods in Hail\n",
    "\n",
    "3) Use Hail to investigate how these methods work on a structured and admixed dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6375472f",
   "metadata": {},
   "source": [
    "# 1. Import and initialize Hail\n",
    "\n",
    "We import Hail, initialize it, and import some plotting tools as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad042803",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ['PYSPARK_SUBMIT_ARGS'] = '--driver-memory 6G pyspark-shell'\n",
    "\n",
    "import hail as hl\n",
    "hl.init()\n",
    "from hail.plot import output_notebook, show\n",
    "from bokeh.models import Slope\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc6b7b7d",
   "metadata": {},
   "source": [
    "# 2. Read HGDP data, add population information, and run PCA\n",
    "\n",
    "This bit is a sped-up review from the first notebook. We want to have the principal component plot as reference for relatedness exploration below!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a08f9b31",
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = hl.read_matrix_table('resources/hgdp.mt')\n",
    "sd = hl.import_table('resources/HGDP_sample_data.tsv',\n",
    "                     key='sample_id',\n",
    "                     impute=True)\n",
    "mt = mt.annotate_cols(sample_idx = hl.int(hl.scan.count()), sample_data = sd[mt.s])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92a5d8cb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "_, scores, _ = hl.hwe_normalized_pca(mt.GT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c7b4243",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(hl.plot.scatter(scores.scores[0], scores.scores[1], \n",
    "                     label=mt.index_cols(scores.key).sample_data.continental_pop,\n",
    "                     size=8,\n",
    "                     xlabel='PC1',\n",
    "                     ylabel='PC2'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "feeea5b0",
   "metadata": {},
   "source": [
    "# 3. Simulate 3 generations of random mating.\n",
    "\n",
    "Hail contains a simple random mating simulator that can be a great tool for exploring relatedness and admixture in populations.\n",
    "\n",
    "The starting dataset is highly **structured** but less **admixed**. We can see distinct ancestral clusters with spaces  between. If we were to make the same plot of a random slice of the human population, the clusters would be connected by individuals falling on the clines between.\n",
    "\n",
    "Here we simulate three rounds of random mating, where each generation creates a number of pairs equal to half its generation size (with replacement, but no self-pairs), and each pair has two children. Each pair is within a single generation; pairs do not span generations. There are no restrictions on pair selection other than that a sample cannot pair with itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "116866b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mt1 = hl.simulate_random_mating(mt,\n",
    "                                pairs_per_generation_multiplier=0.5,\n",
    "                                children_per_pair=2,\n",
    "                                n_rounds=3,\n",
    "                                seed=5).key_cols_by('sample_idx')\n",
    "mt1 = mt1.annotate_cols(**mt.key_cols_by('sample_idx').index_cols(mt1.sample_idx).sample_data).persist()\n",
    "mt1.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a113f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "_, scores1, _ = hl.hwe_normalized_pca(mt1.GT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bb7f66b",
   "metadata": {},
   "source": [
    "#### Take a few moments to explore and discuss the PC plot of the new dataset, and where the new samples fall. What structure is present?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c13de9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined_data = mt1.index_cols(scores1.key)\n",
    "show(hl.plot.scatter(scores1.scores[0], scores1.scores[1],\n",
    "                     label=joined_data.continental_pop,\n",
    "                     hover_fields={'sample_id': joined_data.s},\n",
    "                     size=6,\n",
    "                     xlabel='PC1',\n",
    "                     ylabel='PC2'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50ddec00",
   "metadata": {},
   "source": [
    "# 4. Relatedness estimation\n",
    "\n",
    "Estimation of relatedness between individuals is a core tool in statistical genetics. It is commonly used to verify reported relatedness (both false negatives and false positives), and is an important input to statistical analysis.\n",
    "\n",
    "Some statistical analyses (for instance, GWAS using a simple linear regression) are built around the assumption that variance in genotypes and phenotypes are independent and identically distributed, and an estimated relatedness graph can be used to filter samples to an unrelated set before running these methods.\n",
    "\n",
    "Other methods, like SAIGE, explicitly account for sample covariance to produce statistically sound effect estimation for datasets with related samples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69a3cc38",
   "metadata": {},
   "source": [
    "## Coefficient of kinship versus coefficient of relationship.\n",
    "\n",
    "The kinship coefficient estimated in methods below is defined as the probability that two homologous alleles drawn from each of two individuals are identical by descent. The similar \"coefficient of relationship\", defined as the fraction of genetic material shared identically-by-descent, is equal to twice the kinship coefficient.\n",
    "\n",
    "| Relationship | Kinship Coefficient ($\\phi$) | Coefficient of Relationship ($r$) |\n",
    "| :--- | --- | --- |\n",
    "| Self | 1/2 = 0.5 | 1.0 |\n",
    "| Parent-Child | ??? | ??? |\n",
    "| Full Sibling | ??? | ??? |\n",
    "| Grandparent-Grandchild | 1/8 = 0.125 | 0.25 |\n",
    "| Avuncular Pair | ??? | ??? |\n",
    "| First Cousin | ??? | ??? |\n",
    "\n",
    "### Exercise 1\n",
    "\n",
    "Fill in the missing entries in the above table. Double clicking on a markdown cell, like this one, allows you to edit it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "781182aa",
   "metadata": {},
   "source": [
    "# 5. Identity By Descent\n",
    "\n",
    "All of the models we explore today are so-called \"SNP-based\" methods. These methods operate in two phases. In phase one, each variant is independently used as an indicator of the relatedness of a pair of samples. In phase two, the kinship of a pair of samples is defined as a weighted sum of these independent tests. The three methods we explore today: IBD, KING, and PC-Relate, differ only in their choice of weights!\n",
    "\n",
    "The PLINK or \"identity-by-descent\" estimator of kinship is defined in Purcell, et al. 2007. The definition is rather complex so we do not repeat it here. Instead, consider an equivalent formulation presented Manichaikul, et al. 2010:\n",
    "\n",
    "Let\n",
    "- $X_{im}$ be the number of alternate alleles in the $i$-th sample's genotype at the $m$-th biallelic SNP.\n",
    "- $p_s$ be the population alternate allele frequency of the $m$-th biallelic SNP.\n",
    "- $M$ is the number of biallelic SNPs.\n",
    "\n",
    "The estimator for the kinship between sample $i$ and sample $j$ is given by:\n",
    "\n",
    "$$\n",
    "\\widehat{\\phi_{ij}} = \\frac{1}{2} - \\frac{\\sum_{m} (X_{im} - X_{jm})^2}{\\sum_{m} 8 p_{m} (1 - p_{m})}\n",
    "$$\n",
    "\n",
    "Intuitively: we are counting down from an expectation of monozygotic twins ($\\frac{1}{2}$). If two individuals share all their genotypes, the numerator of the second term goes to zero and $\\widehat{\\phi_{ij}}$ goes to $\\frac{1}{2}$. If two individuals have opposing homozygous genotypes at every site, the numerator goes to $4M$ and the denominator goes to $8M$ (note that the variance tends one if every site tends to an even split between homozygous alternate and homozygous reference)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0be2096",
   "metadata": {},
   "source": [
    "Below we use the relatedness graph produced by `simulate_random_mating` to look up the relatedness for a given pair, and adds a little bit of noise to make the plots look more realistic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd3e865",
   "metadata": {},
   "outputs": [],
   "source": [
    "rel = mt1.index_globals().relatedness\n",
    "def get_rel(s1, s2):\n",
    "    return (hl.case()\n",
    "            .when(s1 == s2, 0.5)\n",
    "            .when(s1 > s2, hl.coalesce(rel.get(s1).get(s2), 0.0))\n",
    "            .default(hl.coalesce(rel.get(s2).get(s1), 0.0))) + hl.rand_norm(0, 0.005) # jitter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bcb757a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ibd = hl.identity_by_descent(mt1.key_cols_by(sample_idx=hl.str(mt1.sample_idx)), min=0, max=1)\n",
    "ibd = ibd.key_by(i=hl.int(ibd.i), j=hl.int(ibd.j))\n",
    "fig = hl.plot.scatter(get_rel(ibd.i, ibd.j) * 2,\n",
    "                      ibd.ibd.PI_HAT, xlabel='True coefficient of relatedness', ylabel='IBD PI_HAT',\n",
    "                      hover_fields={'id1': ibd.i, 'id2': ibd.j})\n",
    "fig.add_layout(Slope(gradient=1, y_intercept=0, line_color='red', line_dash='dashed'))\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3755c78",
   "metadata": {},
   "source": [
    "### Exercise 2: discuss -- How does the estimate deviate from true relatedness? Why?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40d2f4a2",
   "metadata": {},
   "source": [
    "# 6. The KING estimator\n",
    "\n",
    "The estimator for kinship given in Section 5 assumes a single homogenous population. Do you know why?\n",
    "\n",
    "The KING estimator for kinship relaxes this assumption to multiple homogenous populations with no admixture. The estimator was defined by Manichaikul et al. 2010. We present here only the \"between-family\" estimator, which is the only one implemented in Hail.\n",
    "\n",
    "Let\n",
    "\n",
    "- $N^{Aa}_i$ be the number of heterozygote genotypes in sample $i$.\n",
    "\n",
    "The KING estimator for between-family kinship is given by:\n",
    "\n",
    "$$\n",
    "\\widehat{\\phi^{\\textrm{between}}_{ij}} = \\frac{1}{2} - \\frac{\\sum_m (X_{im} - X_{jm})^2}{4\\textrm{min}(N^{Aa}_i, N^{Aa}_j)}\n",
    "$$\n",
    "\n",
    "Intuitively: we are again counting down from an assumption of monozygotic twins. Critically, instead of a sum of the binomial variance for two draws (one for each allele), $2p_m(p_m - 1)$, Manichaikul et al. propose using the minimum count of heterozygous variants between the two samples. Manichaikul et al. justify this proposal by noting that the fraction of heterozygotic genotypes is equal to the expectation of twice the binomial variance of alternate allele frequency:\n",
    "\n",
    "$$\n",
    "N^{Aa}_i / M = \\mathbb{E}_m( 2p_m(1 - p_m) )\n",
    "$$\n",
    "\n",
    "And they justify using the minimum of the two samples' heterozygote counts thusly:\n",
    "\n",
    "> In order to guard against potential estimation inflation due to departure from individual-level HWE, we consider the smaller of the observed heterozygosity rates [...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75d978d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "king = hl.king(mt1.GT).entries()\n",
    "king = king.filter(king.sample_idx != king.sample_idx_1) # remove self-comparisons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2dc9c1e",
   "metadata": {},
   "source": [
    "Now we'll make a plot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd96c6ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = hl.plot.scatter(get_rel(king.sample_idx, king.sample_idx_1),\n",
    "                      king.phi, \n",
    "                      xlabel='True kinship', \n",
    "                      ylabel='KING Phi',\n",
    "                      hover_fields={'id1': king.sample_idx, 'id2': king.sample_idx_1})\n",
    "\n",
    "fig.add_layout(Slope(gradient=1, y_intercept=0, line_color='red', line_dash='dashed'))\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c88ee8ef",
   "metadata": {},
   "source": [
    "### Exercise 3: discuss -- How does the estimate deviate from true relatedness? Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50cd3ad1",
   "metadata": {},
   "source": [
    "# 7. The PC-Relate Estimator\n",
    "\n",
    "When an individual has admixed ancestry, the KING model breaks down. Alleles that are rare in all but one ancestry group are treated as very informative to relatedness. However, these alleles are simply markers of the ancestry group. The PC-Relate method corrects for this.\n",
    "\n",
    "PC-Relate slightly modifies the usual estimator for relatedness: occurrences of population allele frequency are replaced with an \"individual-specific allele frequency\". This modification allows the method to correctly weight an allele according to an individual’s unique ancestry profile.\n",
    "\n",
    "\n",
    "\n",
    "Let \n",
    "- $\\widehat{\\mu_{im}}$ be the \"individual-specific allele frequency\" at a given genetic locus. Define it as one-half the best linear estimate of an individual's genotype using only the provided principal components.\n",
    "- $\\widehat{\\sigma^2_{im}} = \\widehat{\\mu_{im}}(1 - \\widehat{\\mu_{im}})$\n",
    "- $\\widehat{\\sigma_{im}} = \\sqrt{\\widehat{\\sigma^2_{im}}}$\n",
    "\n",
    "The PC-Relate estimator of kinship is then defined as:\n",
    "\n",
    "$$\n",
    "\\widehat{\\phi_{ij}} = \\frac{\\sum_m (g_{im} - 2\\mu_{im}) (g_{jm} - 2\\mu_{jm})}{4 \\sum_m \\widehat{\\sigma_{im}} \\widehat{\\sigma_{jm}}}\n",
    "$$\n",
    "\n",
    "Intuitively: calculate the covariance of the residual genotypes (after \"individual-specific mean\" centering) normalized by the sum of their \"individual-specific allele frequency variances\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3352064b",
   "metadata": {},
   "outputs": [],
   "source": [
    "pcrel = hl.pc_relate(mt1.GT, 0.01, k=4, statistics='kin')\n",
    "pcrel = pcrel.key_by(i=pcrel.i.sample_idx, j=pcrel.j.sample_idx)\n",
    "fig = hl.plot.scatter(get_rel(pcrel.i, pcrel.j),\n",
    "                      pcrel.kin, xlabel='True kinship', ylabel='PC-Relate Kinship',\n",
    "                      hover_fields={'id1': pcrel.i, 'id2': pcrel.j})\n",
    "fig.add_layout(Slope(gradient=1, y_intercept=0, line_color='red', line_dash='dashed'))\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cc8362",
   "metadata": {},
   "source": [
    "### Exercise 4: detective work. Investigate the relationship between the individuals with kinship coefficient ~0.375.\n",
    "\n",
    "Some useful information is encoded in the column fields of `mt`. You can show specific samples by editing the below code.\n",
    "    \n",
    "If you've got time after finishing this, do the same for a pair in the cluster with relatedness ~0.185!\n",
    "\n",
    "To get you started, here's an example of how to interrogate three random samples, sample 50, 100, and 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "407b526c",
   "metadata": {},
   "outputs": [],
   "source": [
    "mt1.filter_cols(hl.literal([50, 100, 1000]).contains(mt1.sample_idx)).cols().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7e41b9b",
   "metadata": {},
   "source": [
    "### Exercise 5: model interrogation. \n",
    "\n",
    "PC-Relate uses an explicit `k` term to control how many principal components are used to in individual allele frequency predictions. Your job is to copy or edit the code cell starting with `pcrel = `, and run it with different `k` values to interrogate how the number of principal components included affects the results.\n",
    "\n",
    "What **k** seems best? What happens when **k** is small? Large?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7688c2d5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}