Genotype_dimred_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimension reduction of genotype data\n",
    "\n",
    "This notebook is intended to explore population structure through dimension reduction of genotype data from the Thousand Genomes Project (1KGP). Data is available here:\n",
    "\n",
    "* ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/hd_genotype_chip/\n",
    "* ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/\n",
    "\n",
    "We will use the following files:\n",
    "* `ALL.wgs.nhgri_coriell_affy_6.20140825.genotypes_has_ped.vcf.gz`\n",
    "* `affy_samples.20141118.panel`\n",
    "* `20131219.populations.tsv`\n",
    "\n",
    "t-SNE is available in sklearn but can be quite slow. Multi-core t-SNE and UMAP are available github:\n",
    "* https://github.com/lmcinnes/umap\n",
    "* https://github.com/DmitryUlyanov/Multicore-TSNE/\n",
    "\n",
    "To install UMAP, run one of the following:\n",
    "\n",
    "```conda install -c conda-forge umap-learn```\n",
    "\n",
    "```pip install umap-learn```\n",
    "\n",
    "If neither works, please visit the UMAP github page for more detailed directions.\n",
    "\n",
    "Papers for t-SNE and UMAP can be found, respectively, at the following links:\n",
    "* https://lvdmaaten.github.io/publications/papers/JMLR_2008.pdf\n",
    "* https://arxiv.org/pdf/1802.03426.pdf\n",
    "\n",
    "This code was written by Alex Diaz-Papkovich and Simon Gravel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries. \n",
    "# Generate images in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import collections\n",
    "from collections import defaultdict\n",
    "import gzip\n",
    "import itertools\n",
    "import numpy as np\n",
    "import os\n",
    "import time\n",
    "\n",
    "from ipywidgets import interact\n",
    "import bokeh\n",
    "import bokeh.io\n",
    "from bokeh.io import push_notebook\n",
    "from bokeh.plotting import figure, show, save, output_notebook, output_file\n",
    "\n",
    "# Import colour palettes for later on\n",
    "from bokeh.palettes import Category20b\n",
    "from bokeh.palettes import Purples\n",
    "from bokeh.palettes import Greens\n",
    "from bokeh.palettes import YlOrBr\n",
    "from bokeh.palettes import YlOrRd\n",
    "from bokeh.palettes import PuOr\n",
    "from bokeh.palettes import RdGy\n",
    "\n",
    "# Dimension reduction tools\n",
    "from sklearn.decomposition import PCA as PCA\n",
    "from sklearn.manifold import TSNE\n",
    "import umap "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify your parent directory. This is where your 1KGP files are stored\n",
    "data_dir = '@@@@'\n",
    "\n",
    "# These are the names of the files we use\n",
    "vcf_name = 'ALL.wgs.nhgri_coriell_affy_6.20140825.genotypes_has_ped.vcf.gz'\n",
    "pop_desc_name = '20131219.populations.tsv'\n",
    "pop_file_name = 'affy_samples.20141118.panel'\n",
    "\n",
    "vcf_file = os.path.join(data_dir, vcf_name)\n",
    "population_description_file = os.path.join(data_dir, pop_desc_name)\n",
    "population_file = os.path.join(data_dir, pop_file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to read in the SNP data. Assign every SNP a value of {0,1,2} relative to reference genome.\n",
    "from collections import Counter\n",
    "\n",
    "class snp(object):\n",
    "\n",
    "    def __init__(self, line, select=False, autosome_only =True):\n",
    "        \"\"\"The initialization method takes in a line from the vcf file, as a string, \n",
    "        and records the relevant information. \n",
    "        line: a string from a vcf file\n",
    "        select: a list of positions of individuals to be analyzed, where positions run from 0 to \n",
    "        nInd-1, the number of individuals\n",
    "        \"\"\" \n",
    "        \n",
    "        split_line = line.split()  #  First break down the line into a list of each field\n",
    "        \n",
    "        self.failed = False  # A label that we will set to True if something goes wrong.\n",
    "        \n",
    "        if line.startswith('#'):\n",
    "            self.failed = True\n",
    "            self.failure_cause = \"line was a header line, not a snp\"\n",
    "            return\n",
    "        \n",
    "        if len(split_line)<=5:\n",
    "            self.failed = True\n",
    "            self.failure_cause = \"incorrectly formatted line, should have at least 5 fields \" + line\n",
    "            return\n",
    "          \n",
    "        self.chrom = split_line[0]\n",
    "        if autosome_only:\n",
    "            if self.chrom not in [\"%d\" % (i,) for i in range(1,23)]:\n",
    "                self.failed = True\n",
    "                self.failure_cause = \"not recognized as an autosome while autosome_only set to True\"\n",
    "                return\n",
    "        \n",
    "        self.chrom = int(split_line[0]) # Chromosome (numbered)\n",
    "        self.position = int(split_line[1])  # The coordinates of the snp\n",
    "        self.rid = split_line[2] # Name/Record ID\n",
    "        self.ref_allele = split_line[3]\n",
    "        self.alt_allele = split_line[4] # The alterate allele according to the vcf; also a string \n",
    "        # Only accept snps in ACGT. \n",
    "        if self.ref_allele not in [\"A\",\"C\",\"G\",\"T\"] or self.alt_allele not in [\"A\",\"C\",\"G\",\"T\"]:\n",
    "            self.failed = True\n",
    "            self.failure_cause = \"ref or alt not in ACGT\"\n",
    "            return\n",
    "        self.filter = split_line[6]  # See vcf format specifications for the interpretation of \n",
    "                                    # the filter field\n",
    "        if self.filter not in ['PASS', '.'] :  # PASS indicates a SNP that passed all QC filters.\n",
    "            self.failed = True\n",
    "            self.failure_cause = self.filter\n",
    "            return\n",
    "              \n",
    "        self.genotype_strings = split_line[9:]\n",
    "\n",
    "        # Prepare a list that will contain the transformed genotypes. \n",
    "        # Since we already know how long the list will be, it makes sense \n",
    "        # to create an array of zeros of the same length as self.gtypes, \n",
    "        \n",
    "        self.genotype_array = np.zeros(len(self.genotype_strings), dtype = np.int8)             \n",
    "\n",
    "        # Count the number of each genotype. \n",
    "        # There may be different strings giving the same genotype so we increment the \n",
    "        # counts found so far for the genotype by the number of times the  \n",
    "        # For example, \"0/0\" and \"0\\0\" give homref, and \"0|1\" and \"1|0\" give het\n",
    "        \n",
    "        n_missing = 0\n",
    "        for index,genotype_string in enumerate(self.genotype_strings):\n",
    "            if genotype_string == './.':\n",
    "                n_missing +=1 \n",
    "                self.genotype_array[index]=-1\n",
    "                continue # missing data will be left as 0\n",
    "            allele_0 = genotype_string[0] # Get the first allele (as a string)\n",
    "            allele_1 = genotype_string[2]\n",
    "            if (allele_0=='1' and allele_1=='1'): # Use rstrip because windows machines will occasionally have extra \\n\n",
    "                self.genotype_array[index]=2\n",
    "            elif ((allele_0=='0' and allele_1=='1') or (allele_0=='1' and allele_1=='0')):\n",
    "                self.genotype_array[index]=1   \n",
    "            elif (allele_0=='0' and allele_1=='0'):\n",
    "                # The array was initialized to zero, so nothing to do here!\n",
    "                continue\n",
    "            else:\n",
    "                print((\"unknown genotype\", genotype_string))\n",
    "                self.failed=True\n",
    "                self.failedreason=\"unknown genotype\"\n",
    "                return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following step imports the genotype data. It is not particularly efficient so will take a few minutes even if we skip some of the lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the number of lines to skip to avoid storing every line in memory\n",
    "number_of_lines_to_skip = 1000\n",
    "\n",
    "start_time = time.time()\n",
    "\n",
    "genotype_matrix = []  # Will contain our numerical genotype matrix. \n",
    "genotype_positions = []\n",
    "genotype_names = []\n",
    "x = 0\n",
    "error_count = 0\n",
    "\n",
    "with gzip.open(vcf_file,'rt') as f:\n",
    "    count = 0\n",
    "    for line in f:\n",
    "        count+=1\n",
    "        if count % number_of_lines_to_skip == 0:\n",
    "            if line.startswith(\"#\") or snp(line).failed:\n",
    "                if snp(line).failure_cause != \"line was a header line, not a snp\":\n",
    "                    error_count += 1\n",
    "                    if x < 10:\n",
    "                        print('Failed: ' + snp(line).failure_cause)\n",
    "                        x+=1\n",
    "                continue\n",
    "            \n",
    "            return_snp = snp(line)\n",
    "            genotype_matrix.append(return_snp.genotype_array)\n",
    "            genotype_names.append(return_snp.rid)\n",
    "            genotype_positions.append([return_snp.chrom, return_snp.position])\n",
    "\n",
    "end_time = time.time()\n",
    "            \n",
    "print(\"Run time in seconds: \" + str(end_time - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Transpose the matrix\n",
    "transposed_genotype_matrix = np.array(genotype_matrix).transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code imports auxiliary data (populations, continent, descriptive data, colouring, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "population_by_individual = defaultdict(int)\n",
    "individuals_by_population = defaultdict(list)  # A dictionary containing all the individuals in a given population\n",
    "\n",
    "for line in open(population_file,'r'):\n",
    "    split_line = line.split()\n",
    "    if split_line[0] == 'sample':  # header line\n",
    "        continue\n",
    "\n",
    "    sample_name = split_line[0]\n",
    "    population_name = split_line[1]\n",
    "    population_by_individual[sample_name] = population_name\n",
    "    individuals_by_population[population_name].append(sample_name) \n",
    "\n",
    "populations = list(individuals_by_population.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The path to the place where you put the population name file.\n",
    "name_by_code = {}  # A dictionary giving the full name of each population code\n",
    "pop_by_continent = {}  # A dictionary giving the code of each population within a continent  \n",
    "continent_by_population = {}  # A dictionary giving the continent for each population code\n",
    "for line in open(population_description_file,'r'):\n",
    "    split_line = line.split('\\t')\n",
    "    if split_line[0] in ['Population Description','Total','']:  # header or footer\n",
    "        continue\n",
    "    name_by_code[split_line[1]] = split_line[0]\n",
    "    continent_by_population[split_line[1]] = split_line[2]\n",
    "    try: \n",
    "        pop_by_continent[split_line[2]].append(split_line[1])\n",
    "    except KeyError:\n",
    "        pop_by_continent[split_line[2]] = [split_line[1]]\n",
    "\n",
    "continents = list(pop_by_continent.keys()) \n",
    "    \n",
    "    \n",
    "# Populations listed by continent\n",
    "pops=[]\n",
    "for continent in continents:\n",
    "    pops.extend(pop_by_continent[continent])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign colours to each population, roughly themed according to continent\n",
    "# The Category20b palette has a bunch of groups of 4 shades in the same colour range\n",
    "color_dict = {}\n",
    "for i, cont in enumerate(continents): \n",
    "    for j, pop in enumerate(pop_by_continent[cont]):\n",
    "        color_dict[pop] = Category20b[20][4*i+j%4]\n",
    "\n",
    "# Colour palette above only really supports groups of 4 so we have to manually specify a few colours for the 5th/6th\n",
    "# members of a group\n",
    "\n",
    "color_dict['CHS'] = Purples[9][4]# purple\n",
    "color_dict['STU'] = Greens[9][6] # green\n",
    "color_dict['LWK'] = PuOr[11][-1] # brown\n",
    "color_dict['MSL'] = PuOr[11][-2] # rusty brown\n",
    "color_dict['YRI'] = PuOr[11][-3] # cappucino w/ extra milk (stirred)\n",
    "color_dict['CEU'] = RdGy[11][-3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for line in gzip.open(vcf_file,'rt'):\n",
    "    if line.startswith(\"#\"):\n",
    "        if not line.startswith(\"##\"):\n",
    "            # Extract the individuals for the population, as a list of strings\n",
    "            # Windows users may have trailing \\n characters\n",
    "            individuals = line.split()[9:]\n",
    "            # Once we've extracted the individuals, we can exit the loops. \n",
    "            break\n",
    "\n",
    "# Build a list of populations for each indiviudal in the vcf file\n",
    "lspop = []\n",
    "for ind in individuals:\n",
    "    pop = population_by_individual[ind]\n",
    "    if pop == 0:\n",
    "        lspop.append(\"missing\")\n",
    "    else:\n",
    "        lspop.append(pop)\n",
    "\n",
    "        \n",
    "indices_of_population_members = defaultdict(list)\n",
    "\n",
    "for index,individual in enumerate(individuals):\n",
    "    try:\n",
    "        indices_of_population_members[population_by_individual[individual]].append(index)\n",
    "    except KeyError: # We do not have population info for this individual\n",
    "        continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count = 0\n",
    "\n",
    "for p in pop_by_continent:\n",
    "    count+=len(pop_by_continent[p])\n",
    "    \n",
    "print(count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imax = 0\n",
    "imin = 200\n",
    "for i in indices_of_population_members:\n",
    "    imax = max(len(indices_of_population_members[i]),imax)\n",
    "    imin = min(len(indices_of_population_members[i]),imin)\n",
    "    \n",
    "print(imax, imin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimension reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we generate our prinicipal component projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the PC axes\n",
    "# This will take a few minutes\n",
    "pca_full = PCA().fit(transposed_genotype_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project onto the PC axes\n",
    "proj_pca = pca_full.fit_transform(transposed_genotype_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code uses non-PCA methods. We explore t-SNE and UMAP specifically. For t-SNE you can work in 2 or 3 dimensions. You can also specify a variety of hyperparameters (e.g. perplexity to control \"tightness\" of clusters) but the defaults are sufficient for our purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project the genotype data matrix to two dimensions via t-SNE. This may take several minutes to run.\n",
    "proj_tsne_gt = TSNE(n_components = 2).fit_transform(transposed_genotype_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "UMAP can project to an arbitrary number of dimensions, but for visualizations we stick to 2. Three-dimensional projections can be interesting too though! The developers recommend using parameter values of 5 to 50 neighbours and a minimum distance between 0.001 and 0.5. We use the default number of neighbours (15) and a minimum distance of 0.5 for clarity in visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project the genotype data matrix to two dimensions via UMAP\n",
    "proj_umap_gt = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.5).fit_transform(transposed_genotype_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what happens when we project the top principal components rather than the genotype data directly? Using more components adds more data, but can also hit diminishing returns and add noise. Feel free to tweak the number of components from as low as 2 to as high as you'd like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of principal components to use\n",
    "n_pc = 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project the principal components via t-SNE to 2 dimensions.\n",
    "proj_tsne_pca = TSNE(n_components=2).fit_transform(proj_pca[:,:n_pc])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project the principal components via UMAP to 2 dimensions.\n",
    "proj_umap_pca = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.5).fit_transform(proj_pca[:,:n_pc])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pc = 25\n",
    "proj_umap_pca_2 = umap.UMAP(n_components=2, n_neighbors=15, min_dist=0.5).fit_transform(proj_pca[:,:n_pc])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate interactive HTML files\n",
    "\n",
    "In this section we will generate interactive HTML files that let you select or deselect populations in the projections. You can use the `save()` command to create a file without viewing it or the `show()` command to create and open the file. By default, the files will be created in the directory of this notebook.\n",
    "\n",
    "The three variables to select are the projection to use and the two dimensions you would like to plot against each other. By default we use dimensions 1 and 2 (indexed as 0 and 1 in Python). In PCA we have many dimensions to work with but it can be interesting to look at, say, dimensions 4 vs 1 or 3 vs 2. Our t-SNE and UMAP projections are only 2D so if you want to go into higher dimensions you'll have to specify that and re-run the code. t-SNE only projects to at most 3 dimensions. UMAP projects to an arbitrary number of dimensions.\n",
    "\n",
    "Note: The HTML files will be 800 high by 1350 wide - you may need to change this to fit into your screen. Otherwise, you can use your browser's native zoom/view options to make things fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the PCA projection\n",
    "dset = proj_pca\n",
    "\n",
    "# Select the dimensions in the projection\n",
    "dim1 = 0\n",
    "dim2 = 1\n",
    "\n",
    "p = figure(plot_width=1350, plot_height=800)\n",
    "p.title.text = 'PCA projection: PC' + str(dim1+1) + ' vs PC' + str(dim2+1)\n",
    "\n",
    "for cont in continents: \n",
    "    for pop in pop_by_continent[cont]:\n",
    "        projections_within_population = dset[indices_of_population_members[pop]]\n",
    "        p.circle(projections_within_population[:,dim1], projections_within_population[:,dim2], \n",
    "                 legend=name_by_code[pop], color = color_dict[pop])\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "output_file(\"interactive_pca.html\", title=\"PCA projection\")\n",
    "\n",
    "#save(p)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the UMAP projection of the genotype data\n",
    "dset = proj_umap_gt\n",
    "\n",
    "dim1 = 0\n",
    "dim2 = 1\n",
    "\n",
    "p = figure(plot_width=1350, plot_height=800)\n",
    "p.title.text = 'UMAP projection of genotype data'\n",
    "\n",
    "for cont in continents: \n",
    "    for pop in pop_by_continent[cont]:\n",
    "        projections_within_population = dset[indices_of_population_members[pop]]\n",
    "        p.circle(projections_within_population[:,dim1], projections_within_population[:,dim2], \n",
    "                 legend=name_by_code[pop], color = color_dict[pop])\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "output_file(\"interactive_umap.html\", title=\"UMAP projection\")\n",
    "\n",
    "#save(p)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the tSNE projection of the genotype data\n",
    "dset = proj_tsne_gt\n",
    "\n",
    "dim1 = 0\n",
    "dim2 = 1\n",
    "\n",
    "p = figure(plot_width=1350, plot_height=800)\n",
    "p.title.text = 't-SNE projection of genotype data'\n",
    "\n",
    "for cont in continents: \n",
    "    for pop in pop_by_continent[cont]:\n",
    "        projections_within_population = dset[indices_of_population_members[pop]]\n",
    "        p.circle(projections_within_population[:,dim1], projections_within_population[:,dim2], \n",
    "                 legend=name_by_code[pop], color = color_dict[pop])\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "output_file(\"interactive_tsne.html\", title=\"t-SNE projection\")\n",
    "\n",
    "#save(p)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the UMAP projection of the PCs\n",
    "dset = proj_umap_pca\n",
    "\n",
    "dim1 = 0\n",
    "dim2 = 1\n",
    "\n",
    "p = figure(plot_width=1350, plot_height=800)\n",
    "p.title.text = 'UMAP projection of principal components'\n",
    "\n",
    "for cont in continents: \n",
    "    for pop in pop_by_continent[cont]:\n",
    "        projections_within_population = dset[indices_of_population_members[pop]]\n",
    "        p.circle(projections_within_population[:,dim1], projections_within_population[:,dim2], \n",
    "                 legend=name_by_code[pop], color = color_dict[pop])\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "output_file(\"interactive_pca_umap.html\", title=\"PCA-UMAP projection\")\n",
    "\n",
    "#save(p)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the UMAP projection of the PCs\n",
    "dset = proj_umap_pca_2\n",
    "\n",
    "dim1 = 0\n",
    "dim2 = 1\n",
    "\n",
    "p = figure(plot_width=1350, plot_height=800)\n",
    "p.title.text = 'UMAP projection of principal components'\n",
    "\n",
    "for cont in continents: \n",
    "    for pop in pop_by_continent[cont]:\n",
    "        projections_within_population = dset[indices_of_population_members[pop]]\n",
    "        p.circle(projections_within_population[:,dim1], projections_within_population[:,dim2], \n",
    "                 legend=name_by_code[pop], color = color_dict[pop])\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "output_file(\"interactive_pca_umap.html\", title=\"PCA-UMAP projection\")\n",
    "\n",
    "#save(p)\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the tSNE projection of the PCs\n",
    "dset = proj_tsne_pca\n",
    "\n",
    "dim1 = 0\n",
    "dim2 = 1\n",
    "\n",
    "p = figure(plot_width=1350, plot_height=800)\n",
    "p.title.text = 't-SNE projection of principal components'\n",
    "\n",
    "for cont in continents: \n",
    "    for pop in pop_by_continent[cont]:\n",
    "        projections_within_population = dset[indices_of_population_members[pop]]\n",
    "        p.circle(projections_within_population[:,dim1], projections_within_population[:,dim2], \n",
    "                 legend=name_by_code[pop], color = color_dict[pop])\n",
    "\n",
    "p.legend.location = \"top_left\"\n",
    "p.legend.click_policy=\"hide\"\n",
    "\n",
    "output_file(\"interactive_pca_tsne.html\", title=\"PCA-TSNE projection\")\n",
    "\n",
    "#save(p)\n",
    "show(p)"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}