mining_pd_genetics_imaging_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mining genetic, transcriptomic and imaging data in Parkinson’s disease - 1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib_venn import venn3\n",
    "from scipy.stats import norm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import subprocess\n",
    "import math\n",
    "import sys\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Required software\n",
    "\n",
    "- Plink (v1.9b)\n",
    "\n",
    "#### Data description\n",
    "PPMI genoptyping data are stored in PLINK file format. They are stored in two genotyping datasets:\n",
    "\n",
    "- IMMUNO\n",
    "\n",
    "- NEUROX\n",
    "\n",
    "IMMUNO dataset targets genetic loci that are known to be associated with major autoimmune and neuroinflammatory diseases. It contains 196,524 genetic variants: 718 indels and 195,806 SNPs. Among these genetic variants, 1920 SNPs replicates PD associated genetic loci.\n",
    "\n",
    "NeuroX dataset targets ~240,000 exonic genomic variants and ~24,000 custom variants, involved in neurological diseases.\n",
    "\n",
    "Both datasets provides 3 files:\n",
    "- BED\n",
    "- BIM\n",
    "- FAM\n",
    "\n",
    "#### PLINK BIM\n",
    "Plink BIM is a text file with no header, and one line per variant with the following six fields:\n",
    "- Chromosome code (either an integer, or 'X'/'Y'/'XY'/'MT'; '0' indicates unknown) or name\n",
    "- Variant identifier\n",
    "- Position in morgans or centimorgans\n",
    "- Base-pair coordinate (1-based)\n",
    "- Allele 1 (usually minor)\n",
    "- Allele 2 (usually major)\n",
    "\n",
    "#### PLINK FAM\n",
    "Plink FAM is a text file with no header, and one line per sample with the following six fields:\n",
    "- Family ID ('FID')\n",
    "- Within-family ID ('IID'; cannot be '0')\n",
    "- Within-family ID of father ('0' if father isn't in dataset)\n",
    "- Within-family ID of mother ('0' if mother isn't in dataset)\n",
    "- Sex code ('1' = male, '2' = female, '0' = unknown)\n",
    "- Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)\n",
    "\n",
    "#### PLINK BED\n",
    "Plink BED file is the primary representation of genotype calls at biallelic variants. It must be accompanied by .bim and .fam files. Plink BED file is binary, its human-readable version is the PED file (for further information see [here](https://www.cog-genomics.org/plink/1.9/formats#ped))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genotyping_path = \"data/genotyping/\"\n",
    "pheno_path = \"data/pheno/\"\n",
    "immuno_path = os.path.join(genotyping_path, \"IMMUNO\")\n",
    "neurox_path = os.path.join(genotyping_path, \"NEUROX\")\n",
    "hapmap_path = os.path.join(genotyping_path, \"HapMap3_b37\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A first glimpse at data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "immuno_bim = pd.read_csv(\n",
    "    os.path.join(immuno_path, \"IMMUNO.bim\"), \n",
    "    header=None,\n",
    "    sep=\"\\t\"\n",
    ")\n",
    "immuno_bim.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neurox_bim = pd.read_csv(\n",
    "    os.path.join(neurox_path, \"NEUROX.bim\"),\n",
    "    header=None,\n",
    "    sep=\"\\t\"\n",
    ")\n",
    "neurox_bim.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IMMUNO variant rs9729550 should be mapped at position 1135242 in hg19 release. <br>\n",
    "If it was mapped on hg18 release it should be mapped at position 1125105. <br>\n",
    "Therefore, IMMUNO SNPs have been mapped on hg18 assembly!\n",
    "\n",
    "Moreover, NEUROX SNPs have been mapped on hg19 release.\n",
    "\n",
    "This means that we need to lift IMMUNO SNPs positions. Let's do it with Plink and UCSC genome browser tools!\n",
    "\n",
    "First, we need to encode IMMUNO's BIM in a format that the UCSC genome browser can understand: UCSC BED format: \n",
    "\n",
    "```chrom start stop name```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# minimal UCSC BED file requires\n",
    "#\n",
    "# CHROM    START    STOP    FEATURE_ID \n",
    "\n",
    "tmp_chr = immuno_bim.iloc[:,0]  # chromosomes\n",
    "tmp_pos = immuno_bim.iloc[:,3]  # SNP positions\n",
    "tmp_name = immuno_bim.iloc[:,1] # SNP IDs\n",
    "\n",
    "# append 'chr' in front of chromosome numbers (required by liftOver)\n",
    "tmp_chr_str = [''.join([\"chr\", str(c)]) for c in tmp_chr.values]\n",
    "\n",
    "# replace chr23, chr24, chr25 and chr26 with chrX, chrY, chrX and chrMT\n",
    "tmp_chr_str2 = [c.replace('23', 'X') for c in tmp_chr_str]\n",
    "tmp_chr_str3 = [c.replace('24', 'Y') for c in tmp_chr_str2]\n",
    "tmp_chr_str4 = [c.replace('25', 'X') for c in tmp_chr_str3]\n",
    "tmp_chr_str5 = [c.replace('26', 'MT') for c in tmp_chr_str4]\n",
    "\n",
    "# write the resulting UCSC BED file\n",
    "bed = pd.concat(\n",
    "    [pd.DataFrame(tmp_chr_str5), tmp_pos, tmp_pos + 1, tmp_name],\n",
    "    axis=1\n",
    ")\n",
    "print(bed.head(n=5))\n",
    "print(bed.tail(n=5))\n",
    "bed.to_csv(\n",
    "    os.path.join(immuno_path, \"IMMUNO_tolift.bed\"), sep=\"\\t\",\n",
    "    index=False, header=False\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now use [LiftOver](http://genome.ucsc.edu/cgi-bin/hgLiftOver) tool from UCSC genome browser to lift IMMUNO SNPs positions from hg18 to hg19 genome release. Some SNPs will not be converted (four), but we can ignore it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bed_lifted = pd.read_csv(os.path.join(immuno_path, \"IMMUNO_lifted.bed\"), sep=\"\\t\", header=None)\n",
    "\n",
    "# we keep succesfully mapped SNPs\n",
    "good_snps = bed_lifted.iloc[:,3]\n",
    "good_snps.to_csv(\n",
    "    os.path.join(immuno_path, \"good_snps.txt\"),\n",
    "    index=False, header=False\n",
    ")\n",
    "\n",
    "# call plink from command line\n",
    "!plink --bfile {os.path.join(immuno_path, \"IMMUNO\")} --update-map {os.path.join(immuno_path, \"IMMUNO_lifted.bed\")} 2 4 --make-bed --extract {os.path.join(immuno_path, \"good_snps.txt\")} --out {os.path.join(immuno_path, \"IMMUNO_hg19\")}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "immuno_hg19_bim = pd.read_csv(\n",
    "    os.path.join(immuno_path, \"IMMUNO_hg19.bim\"),\n",
    "    sep=\"\\t\",\n",
    "    header=None\n",
    ")\n",
    "immuno_hg19_bim.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now rs9729550 is correctly mapped at position 1135242!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we will create a single dataset combining IMMUNO and NEUROX datasets. Plink will help us doing this!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "immuno_fam = pd.read_csv(\n",
    "    os.path.join(immuno_path, \"IMMUNO_hg19.fam\"), \n",
    "    sep=\" \", \n",
    "    header=None\n",
    ")\n",
    "neurox_fam = pd.read_csv(\n",
    "    os.path.join(neurox_path, \"NEUROX.fam\"), \n",
    "    sep=\" \", \n",
    "    header=None\n",
    ")\n",
    "# recover subjects IDs\n",
    "immuno_subj = immuno_fam.iloc[:,1].tolist()\n",
    "neurox_subj = neurox_fam.iloc[:,1].tolist()\n",
    "\n",
    "# retrieve subjects with data in both IMMUNO and NEUROX \n",
    "common_subj = set(immuno_subj).intersection(set(neurox_subj))\n",
    "\n",
    "# write common subjects to file\n",
    "common_subj_fn = \"data/genotyping/common_subj.txt\"\n",
    "pd.DataFrame(list(zip(common_subj, common_subj)), columns=[\"FID\", \"IID\"]).to_csv(\n",
    "    common_subj_fn,\n",
    "    header=False,\n",
    "    index=False,\n",
    "    sep=\" \"\n",
    ")\n",
    "\n",
    "# call plink on IMMUNO\n",
    "!plink --bfile {os.path.join(immuno_path, \"IMMUNO_hg19\")} --keep {common_subj_fn} --make-bed --out {os.path.join(immuno_path, \"IMMUNO_common\")}\n",
    "# call plink on NEUROX\n",
    "!plink --bfile {os.path.join(neurox_path, \"NEUROX\")} --keep {common_subj_fn} --make-bed --out {os.path.join(neurox_path, \"NEUROX_common\")}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IMMUNO and NEUROX SNPs were not named using standard variant IDs (e.g. RsIDs). Therefore, ther can be identical SNPs called with different IDs in the two datasets. <br>\n",
    "To solve this problem we rename each SNP in the two datasets as:<br>\n",
    "```chrom:pos_A1_A2```\n",
    "\n",
    "For simplicity we provide the new BIM, BED and FAM files renamed (```IMMUNO_renamed``` and ```NEUROX_renamed``` files)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "immuno_bim = pd.read_csv(\n",
    "    os.path.join(immuno_path, \"IMMUNO_renamed.bim\"),\n",
    "    sep=\"\\t\",\n",
    "    header=None\n",
    ")\n",
    "immuno_bim.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the names of SNPs only available in IMMUNO and write \n",
    "# the SNPs appearing uniquely in IMMUNO to a TXT file\n",
    "neurox_bim = pd.read_csv(\n",
    "    os.path.join(neurox_path, \"NEUROX_renamed.bim\"),\n",
    "    sep=\"\\t\",\n",
    "    header=None\n",
    ")\n",
    "only_immuno_snps_fn = \"data/genotyping/onlyIMMUNOsnps.txt\"\n",
    "neurox_snpids_set = set(neurox_bim.iloc[:,1].to_list())\n",
    "immuno_bim[~immuno_bim[1].isin(neurox_snpids_set)].iloc[:,1].to_csv(\n",
    "    only_immuno_snps_fn, header=False, index=False\n",
    ")\n",
    "\n",
    "# call plink\n",
    "!plink --bfile {os.path.join(immuno_path, \"IMMUNO_renamed\")} --extract {only_immuno_snps_fn} --make-bed --out {os.path.join(immuno_path, \"IMMUNO_uniquesnps\")}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now merge the two datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmi_merge_fn = \"data/genotyping/PPMI_merge\"\n",
    "\n",
    "!plink --bfile {os.path.join(neurox_path, \"NEUROX_renamed\")} --bmerge {os.path.join(immuno_path, \"IMMUNO_uniquesnps\")} --make-bed --out {ppmi_merge_fn}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plink produces some warnings on SNPs with different positions and same positions with different SNP names. <br>\n",
    "To make our lives easier, we remove them.\n",
    "\n",
    "Very easy with Unix shell: <br>\n",
    "```\n",
    "cat data/genotyping/PPMI_merge.log | grep 'Warning' | grep 'Variants' | grep -o \"'.*'\" | cut -d' ' -f1,3 | awk '{print $1\"\\n\"$2;}' | tr -d \"'\" | sort -u > data/genotyping/merge_rm_snps.txt\n",
    "```\n",
    "\n",
    "We can now run again Plink excluding the annoying SNPs!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "!plink --bfile {os.path.join(neurox_path, \"NEUROX_renamed\")} --bmerge {os.path.join(immuno_path, \"IMMUNO_uniquesnps\")} --exclude data/genotyping/merge_rm_snps.txt --make-bed --out {ppmi_merge_fn}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our starting dataset is ready, we can move on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality Control"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quality control (QC) is one of the key steps of GWAS analysis. <br>\n",
    "\n",
    "QC ensures that potential errors or artifacts (e.g. contaminations, poor quality of DNA samples, etc.) that could provide wrong or biased results are removed from our data.\n",
    "\n",
    "Let's use Plink to perform QC on PPMI data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Individual and SNP missingness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin by investigating the missingness per individual and per SNPs.\n",
    "\n",
    "We first remove SNPs with missingess rate > 5%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {ppmi_merge_fn} --geno 0.05 --make-bed --out {\"_\".join([ppmi_merge_fn, \"geno\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And remove also those individuals with missingness rate > 5%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"geno\"])} --mind 0.05 --make-bed --out {\"_\".join([ppmi_merge_fn, \"mind\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sex discrepancies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's continue our QC by performing some simple subject level QC. <br>\n",
    "We now remove subjects where provided sex information does not match the genotype inferred sex (could be a data quality issue!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmi_merge_fam = pd.read_csv(\n",
    "    os.path.join(genotyping_path, \"PPMI_merge_mind.fam\"),\n",
    "    sep=\"\\s+\",\n",
    "    header=None\n",
    ")\n",
    "ppmi_merge_fam.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"mind\"])} --check-sex --out {\"_\".join([ppmi_merge_fn, \"sexcheck\"])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmi_checksex = pd.read_csv(\n",
    "    \".\".join([\"_\".join([ppmi_merge_fn, \"sexcheck\"]), \"sexcheck\"]),\n",
    "    sep=\"\\s+\"\n",
    ")\n",
    "ppmi_checksex.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmi_checksex[ppmi_checksex.STATUS==\"PROBLEM\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subj_problem_fn = \"data/genotyping/subjs_toremove.txt\"\n",
    "ppmi_checksex[ppmi_checksex.STATUS == \"PROBLEM\"].iloc[:,0:2].to_csv(\n",
    "    subj_problem_fn,\n",
    "    sep=\" \",\n",
    "    index=False\n",
    ")\n",
    "\n",
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"mind\"])} --remove {subj_problem_fn} --make-bed --out {\"_\".join([ppmi_merge_fn, \"qc_subjs\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, we restrict our analysis to autosomal chromosomes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"qc_subjs\"])} --autosome --make-bed --out {\"_\".join([ppmi_merge_fn, \"autosome\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Minor Allele Frequency (MAF) QC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We continue by removing SNPs with low minor allele frequency. <br>\n",
    "This depends on dataset size, we used < 1% here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"autosome\"])} --maf 0.01 --make-bed --out {\"_\".join([ppmi_merge_fn, \"maf\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hardy-Weinberg equilibrium QC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also remove those SNPs not following Hardy-Weinberg equilibrium (we used a threshold of $1e^{-6}$ here)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"maf\"])} --hwe 1e-6 --make-bed --out {\"_\".join([ppmi_merge_fn, \"hwe\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remove ambigous SNPs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now find A/T and C/G SNPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmi_merge_bim = pd.read_csv(\n",
    "    \"_\".join([ppmi_merge_fn, \"hwe.bim\"]),\n",
    "    sep=\"\\t\",\n",
    "    header=None\n",
    ")\n",
    "bad1 = ppmi_merge_bim[(ppmi_merge_bim.iloc[:,4] == 'A') & (ppmi_merge_bim.iloc[:,5] == 'T')]\n",
    "bad2 = ppmi_merge_bim[(ppmi_merge_bim.iloc[:,4] == 'T') & (ppmi_merge_bim.iloc[:,5] == 'A')]\n",
    "bad3 = ppmi_merge_bim[(ppmi_merge_bim.iloc[:,4] == 'C') & (ppmi_merge_bim.iloc[:,5] == 'G')]\n",
    "bad4 = ppmi_merge_bim[(ppmi_merge_bim.iloc[:,4] == 'G') & (ppmi_merge_bim.iloc[:,5] == 'C')]\n",
    "\n",
    "badSnps = pd.concat([bad1, bad2, bad3, bad4])\n",
    "badSnps_set = set(badSnps.iloc[:,1].values.tolist())\n",
    "\n",
    "badsnps_fn = \"data/genotyping/badsnps.txt\"\n",
    "pd.DataFrame(list(badSnps_set)).to_csv(badsnps_fn, header=False, index=False)\n",
    "\n",
    "# call plink\n",
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"hwe\"])} --exclude {badsnps_fn} --make-bed --out {\"_\".join([ppmi_merge_fn, \"snp\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Genetic ancestry QC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now analyse population genetic stratification.<br>\n",
    "We will merge our dataset with data from the HapMap consortium, who genotyped a lot of people from different populations around the world."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_fn = \"data/genotyping/PPMI_HM\"\n",
    "ppmi_merge_bim = pd.read_csv(\n",
    "    \"_\".join([ppmi_merge_fn, \"snp.bim\"]),\n",
    "    sep=\"\\t\",\n",
    "    header=None\n",
    ")\n",
    "\n",
    "# renaming for HM3 was done as described above\n",
    "hapmap_bim = pd.read_csv(\n",
    "    os.path.join(hapmap_path, \"HM3_b37_renamed.bim\"),\n",
    "    sep=\"\\t\",\n",
    "    header=None\n",
    ")\n",
    "\n",
    "# get SNPs in common between the two datasets\n",
    "A = set(ppmi_merge_bim.iloc[:,1].tolist())\n",
    "B = set(hapmap_bim.iloc[:,1].tolist())\n",
    "hp_snps = A.intersection(B)\n",
    "pd.DataFrame(list(hp_snps)).to_csv(\n",
    "    \"data/genotyping/common_snps.txt\",\n",
    "    index=False,\n",
    "    header=False\n",
    ")\n",
    "\n",
    "# call plink\n",
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"snp\"])} --bmerge {os.path.join(hapmap_path, \"HM3_b37_renamed\")} --extract data/genotyping/common_snps.txt --make-bed --out {hp_fn} "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As during the previous merge we got some warnings from Plink on SNPs from different datasets with different names mappoed at the same position. <br>\n",
    "Let's use again Unix shell to recover them:\n",
    "```\n",
    "cat data/genotyping/PPMI_HM.log | grep 'Warning' | grep 'Variants' | grep -o \"'.*'\" | cut -d' ' -f1,3 | awk '{print $1\"\\n\"$2;}' | tr -d \"'\" | sort -u > data/genotyping/merge_rm_snps.txt\n",
    "```\n",
    "\n",
    "Now we can rerun Plink excluding the annoying SNPs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"snp\"])} --bmerge {os.path.join(hapmap_path, \"HM3_b37_renamed\")} --exclude data/genotyping/merge_rm_snps.txt --extract data/genotyping/common_snps.txt --make-bed --out {hp_fn} "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we remove SNPs with high missingness in the data (```--geno 0.1```) and low minor allele frequency (```--maf 0.05```)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {hp_fn} --geno 0.1 --maf 0.05 --make-bed --out {\"_\".join([hp_fn, \"qc\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use PLINK to compute the first few principal components of the relatedness matrix. <br>\n",
    "This will help us to visualize the population structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([hp_fn, \"qc\"])} --pca 20 --out {\"_\".join([hp_fn, \"qc\"])} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp_pca = pd.read_csv(\n",
    "    \"_\".join([hp_fn, \"qc.eigenvec\"]),\n",
    "    sep=\" \",\n",
    "    header=None\n",
    ")\n",
    "cnames = ['FID','IID']\n",
    "for i in range(1,21):\n",
    "    cnames.append(\"PC\"+str(i))\n",
    "hp_pca.columns = cnames\n",
    "hp_pca.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assign_pop(pop):\n",
    "    if pop.isdigit():\n",
    "        return \"PPMI\"\n",
    "    return pop\n",
    "hp_pca[\"FID\"] = hp_pca.apply(lambda x : assign_pop(x[0]), axis=1)\n",
    "\n",
    "# plot population structure\n",
    "palette = [\n",
    "    \"#000000\", \n",
    "    \"#fbb696\", \n",
    "    \"#913502\", \n",
    "    \"#b48060\", \n",
    "    \"#724a20\", \n",
    "    \"#ffb55f\", \n",
    "    \"#108300\", \n",
    "    \"#006435\", \n",
    "    \"#59dda7\", \n",
    "    \"#508aff\", \n",
    "    \"#53459f\", \n",
    "    \"#cd87ff\"\n",
    "]\n",
    "\n",
    "# PC1 vs PC2\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.scatterplot(data=hp_pca, x=\"PC1\", y=\"PC2\", hue=\"FID\", palette=palette)\n",
    "plt.title(\"Population structure\", size=16)\n",
    "plt.xlabel(\"PC1\", size=14)\n",
    "plt.ylabel(\"PC2\", size=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PC2 vs PC3\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.scatterplot(data=hp_pca, x=\"PC2\", y=\"PC3\", hue=\"FID\", palette=palette)\n",
    "plt.title(\"Population structure\", size=16)\n",
    "plt.xlabel(\"PC2\", size=14)\n",
    "plt.ylabel(\"PC3\", size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now select only those subjects of European (CEU and TSI) ancestry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we measure the distance from the central european (CEU) and italian (TSI) cluster\n",
    "hp_ceu = hp_pca[(hp_pca.FID =='CEU') | (hp_pca.FID =='TSI')]\n",
    "ceu_means = hp_ceu.loc[:,'PC1':'PC5'].apply(np.mean)\n",
    "ceu_sds   = hp_ceu.loc[:,'PC1':'PC5'].apply(np.std)\n",
    "\n",
    "#compute z-scores for PPMI subjects\n",
    "ppmi_pca = hp_pca[hp_pca.FID == \"PPMI\"]\n",
    "ppmi_ceu_z = ppmi_pca.loc[:,'PC1':'PC5'].apply(lambda x: (x - ceu_means)/ceu_sds, axis=1)\n",
    "\n",
    "ppmi_ceu_z.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "keep_ppmi = ppmi_ceu_z.apply(lambda x: abs(x) > 5).apply(np.sum,axis=1) == 0\n",
    "\n",
    "ppmi_color=[\"b\"] * len(keep_ppmi)\n",
    "idx = [i for i, x in enumerate(keep_ppmi.values) if not x]\n",
    "for i in idx:\n",
    "    ppmi_color[i] = \"r\"\n",
    "\n",
    "cuse = hp_pca.apply(lambda x: \"#000000\", axis=1)\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.scatterplot(data=hp_pca, x=\"PC1\", y=\"PC2\", hue=\"FID\", palette=palette)\n",
    "plt.scatter(ppmi_pca.PC1, ppmi_pca.PC2, s=2, c=ppmi_color)\n",
    "plt.title(\"Population structure\", size=16)\n",
    "plt.xlabel(\"PC1\", size=14)\n",
    "plt.ylabel(\"PC2\", size=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmi_pca_ceu = ppmi_pca[keep_ppmi]\n",
    "pd.DataFrame(\n",
    "    zip(ppmi_pca_ceu.IID.tolist(), ppmi_pca_ceu.IID.tolist())\n",
    ").to_csv(\n",
    "    \"data/genotyping/ceu_subjs.txt\",\n",
    "    sep=\" \",\n",
    "    index=False,\n",
    "    header=False\n",
    ")\n",
    "\n",
    "# call plink\n",
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"snp\"])} --keep data/genotyping/ceu_subjs.txt --make-bed --out {\"_\".join([ppmi_merge_fn, \"ceu\"])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now recover the subset of samples with complete genotyping, transcriptomic and neuroimaging data for the following analyses. <br>\n",
    "We also recompute the first 20 PCs of the relatedness matrix which will be used during SNP-phenotype associations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {\"_\".join([ppmi_merge_fn, \"ceu\"])} --keep {os.path.join(pheno_path, \"rnaseq_subjs.txt\")} --make-bed --out {os.path.join(genotyping_path, \"PPMI_merge_final\")} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!plink --bfile {os.path.join(genotyping_path, \"PPMI_merge_final\")}  --pca 20 --out {os.path.join(genotyping_path, \"PPMI_merge_final\")} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}