notebooks/GWAS/GWAS_coat_color.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7a244bb3",
   "metadata": {},
   "source": [
    "# GWAS in the cloud\n",
    "## Overview\n",
    "We adapted the NIH CFDE tutorial from [here](https://training.nih-cfde.org/en/latest/Bioinformatic-Analyses/GWAS-in-the-cloud/background/) and fit it to a notebook. We have greatly simplified the instructions, so if you need or want more details, look at the full tutorial to find out more.\n",
    "\n",
    "Most of this notebook is written in Bash, but expects that you are using a Python kernel, until step 3, plotting where you will need to switch your kernel to R."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3edafe63",
   "metadata": {},
   "source": [
    "## Learning Objectives\n",
    "The goal is to learn how to execute a GWAS analysis in a cloud environment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d7ef396",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "+ You only need access to a Sagemaker notebook environment to run this notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39ee9668",
   "metadata": {},
   "source": [
    "## Get Started"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fbf6304",
   "metadata": {},
   "source": [
    "### Install packages and set up environment\n",
    "\n",
    "#### Download the data\n",
    "use %%bash to denote a bash block. You can also use '!' to denote a single bash command within a Python notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ec900bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir GWAS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c20d2c1-1f0a-441d-9b5e-4ed018ec01ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "curl -LO https://de.cyverse.org/dl/d/E0A502CC-F806-4857-9C3A-BAEAA0CCC694/pruned_coatColor_maf_geno.vcf.gz\n",
    "curl -LO https://de.cyverse.org/dl/d/3B5C1853-C092-488C-8C2F-CE6E8526E96B/coatColor.pheno"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d43ae73",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mv *.gz GWAS\n",
    "mv *.pheno GWAS\n",
    "ls GWAS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37567b92-71f1-4c7d-8cce-c08e82b25557",
   "metadata": {},
   "source": [
    "### Install dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f5032d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# install mamba\n",
    "! curl -L -O https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-$(uname)-$(uname -m).sh\n",
    "! bash Mambaforge-$(uname)-$(uname -m).sh -b -p $HOME/mambaforge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a5bd340",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add to your path\n",
    "import os\n",
    "os.environ[\"PATH\"] += os.pathsep + os.environ[\"HOME\"]+\"/mambaforge/bin\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b219074a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# install everything else\n",
    "! mamba install -y -c bioconda plink vcftools"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3de2fc4c",
   "metadata": {},
   "source": [
    "## Analyze"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "013d960d",
   "metadata": {},
   "source": [
    "### Make map and ped files from the vcf file to feed into plink"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e91c7a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "cd GWAS\n",
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6570875d",
   "metadata": {},
   "outputs": [],
   "source": [
    "! vcftools --gzvcf pruned_coatColor_maf_geno.vcf.gz --plink --out coatColor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9a38761",
   "metadata": {},
   "source": [
    "### Create a list of minor alleles\n",
    "For more info on these terms, look at step 2 at https://training.nih-cfde.org/en/latest/Bioinformatic-Analyses/GWAS-in-the-cloud/analyze/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c868a67",
   "metadata": {},
   "outputs": [],
   "source": [
    "#unzip vcf\n",
    "! vcftools --gzvcf pruned_coatColor_maf_geno.vcf.gz --recode --out pruned_coatColor_maf_geno"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e11f991",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create list of minor alleles\n",
    "! cat pruned_coatColor_maf_geno.recode.vcf | awk 'BEGIN{FS=\"\\t\";OFS=\"\\t\";}/#/{next;}{{if($3==\".\")$3=$1\":\"$2;}print $3,$5;}'  > minor_alleles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cff47e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "! head minor_alleles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56d901c7",
   "metadata": {},
   "source": [
    "### Run quality controls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dafa14a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate missingness per locus\n",
    "! plink --file coatColor --make-pheno coatColor.pheno \"yellow\" --missing --out miss_stat --noweb --dog --reference-allele minor_alleles --allow-no-sex --adjust"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cf5f51b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#take a look at lmiss, which is the per locus rates of missingness\n",
    "! head miss_stat.lmiss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "915bb263",
   "metadata": {},
   "outputs": [],
   "source": [
    "#peek at imiss which is the individual rates of missingness\n",
    "! head miss_stat.imiss"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c11ca71",
   "metadata": {},
   "source": [
    "### Convert to plink binary format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b8f2d7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "! plink --file coatColor --allow-no-sex --dog --make-bed --noweb --out coatColor.binary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e36f6cd7",
   "metadata": {},
   "source": [
    "### Run a simple association step (the GWAS part!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f926ef9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "! plink --bfile coatColor.binary --make-pheno coatColor.pheno \"yellow\" --assoc --reference-allele minor_alleles --allow-no-sex --adjust --dog --noweb --out coatColor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b397d484",
   "metadata": {},
   "source": [
    "### Identify statistical cutoffs\n",
    "This code finds the equivalent of 0.05 and 0.01 p value in the negative-log-transformed p values file. We will use these cutoffs to draw horizontal lines in the Manhattan plot for visualization of haplotypes that cross the 0.05 and 0.01 statistical threshold (i.e. have a statistically significant association with yellow coat color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b94e1e2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "unad_cutoff_sug=$(tail -n+2 coatColor.assoc.adjusted | awk '$10>=0.05' | head -n1 | awk '{print $3}')\n",
    "unad_cutoff_conf=$(tail -n+2 coatColor.assoc.adjusted | awk '$10>=0.01' | head -n1 | awk '{print $3}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f52e97c",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "In this tutorial, plotting is done in R, so at this point you can change your kernel to R in the top right. Wait for it to say 'idle' in the bottom left, then continue. You could also plot using Python native packages and maintain the Python notebook kernel."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "effb5acd",
   "metadata": {},
   "source": [
    "### Install qqman"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60feed89",
   "metadata": {},
   "outputs": [],
   "source": [
    "install.packages('qqman', contriburl=contrib.url('http://cran.r-project.org/'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3f1fcd2",
   "metadata": {},
   "source": [
    "### Run the plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7e8cd2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#make sure you are still CD in GWAS, when you change kernel it may reset to home\n",
    "setwd('GWAS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7946a3a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "require(qqman)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d28ef2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=read.table(\"coatColor.assoc\", header=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e5207be",
   "metadata": {},
   "outputs": [],
   "source": [
    "data=data[!is.na(data$P),]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6330b1e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "manhattan(data, p = \"P\", col = c(\"blue4\", \"orange3\"),\n",
    "          suggestiveline = 12,\n",
    "          genomewideline = 15,\n",
    "          chrlabs = c(1:38, \"X\"), annotateTop=TRUE, cex = 1.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26787d84",
   "metadata": {},
   "source": [
    "In our graph, haplotypes in four parts of the genome (chromosome 2, 5, 28 and X) are found to be associated with an increased occurrence of the yellow coat color phenotype.\n",
    "\n",
    "The top associated mutation is a nonsense SNP in the gene MC1R known to control pigment production. The MC1R allele encoding yellow coat color contains a single base change (from C to T) at the 916th nucleotide."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f6e1ef6",
   "metadata": {},
   "source": [
    "### Conclusion\n",
    "Here we learned how to run a simple GWAS analysis in the cloud"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "044a04d8",
   "metadata": {},
   "source": [
    "## Clean up\n",
    "Make sure you shut down this VM, or delete it if you don't plan to use if further.\n",
    "\n",
    "You can also [delete the buckets](https://docs.aws.amazon.com/AmazonS3/latest/userguide/delete-bucket.html) if you don't want to pay for the data: `aws s3 rb s3://bucket-name --force`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1e7be16",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}