AMIA_Exercise6-PheWAS_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PheWAS Analysis of FTO Variant rs8050136\n",
    "\n",
    "This notebook performs the PheWAS analysis of our extracted variant. The 16 cpu + 104gb ram is sufficient for this single variant analysis, though it will take around 20 minutes per analysis. This PheWAS is the simplest form- we are not stratifying on genomic ancestry here, which may be wise for more scientifically robust anaylses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load libraries and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(PheWAS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(tidyverse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Covariates\n",
    "\n",
    "From exercises 3, 4, and 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covariates = read_csv(\"covariates.csv\")\n",
    "covariates %>% transmute(across(-person_id, \\(x) sample(x))) %>% head()\n",
    "#head(covariates) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variant information\n",
    "\n",
    "From exercises 1 and 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snp = read_delim(\"fto_snp.raw\", delim=\" \")  %>% transmute(person_id=IID, rs8050136_A=`chr16:53782363:C:A_A`)\n",
    "snp %>% transmute(across(-person_id, \\(x) sample(x))) %>% head()\n",
    "#head(snp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Phecode phenotypes\n",
    "\n",
    "From exercise 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phenotypes = read_csv(\"phecode_table.csv\")\n",
    "phenotypes[,1:10] %>% transmute(across(-person_id, \\(x) sample(x))) %>% head()\n",
    "#head(phenotypes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge the tables for analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data = inner_join(inner_join(covariates, snp), phenotypes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the PheWAS\n",
    "\n",
    "This will execute our PheWAS analysis. We will use every phecode phenotype (The names of our phenotypes table, except person_id) and use a few covariates. This analysis uses all 16 cores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = phewas_ext(data, phenotypes = names(phenotypes)[-1], \n",
    "                     genotypes=c(\"rs8050136_A\"),\n",
    "                    covariates=c(\"sex_ab_F\", \"age\", \"pc1\", \"pc2\", \"pc3\"),\n",
    "                    cores=32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Results\n",
    "\n",
    "I find these plot options yield an easy, nice looking plot in the notebook environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=15, repr.plot.height=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phewasManhattan(results, OR.direction=T, title =\"PheWAS against rs8050136_A\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show results tables\n",
    "\n",
    "The phewasDT() function creates a handy interactive table. It won't render in the workbench preview, but is nice when you are looking at results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phewasDT(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Static results table\n",
    "\n",
    "This code gives a simple table view of top results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p.digits=3\n",
    "OR.digits=2\n",
    "results %>% addPhecodeInfo() %>% filter(!is.na(p)) %>% arrange(p) %>% \n",
    "        transmute(phenotype, description, snp, OR = round(OR, \n",
    "            OR.digits), p = signif(p, p.digits), n_total, n_cases, \n",
    "            n_controls) %>% head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BMI-Adjusted PheWAS Analysis\n",
    "\n",
    "We can also adjust for BMI in our analysis. The next section does the same type of analysis, but adjusted for BMI. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_adj = phewas_ext(data, phenotypes = names(phenotypes)[-1], \n",
    "                     genotypes=c(\"rs8050136_A\"),\n",
    "                    covariates=c(\"sex_ab_F\", \"age\", \"pc1\", \"pc2\", \"pc3\", \"BMI\"),\n",
    "                    cores=32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phewasManhattan(results_adj, OR.direction=T, title =\"PheWAS against rs8050136_A, BMI-Adjusted\",\n",
    "               y.axis.interval=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phewasDT(results_adj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_adj %>% addPhecodeInfo() %>% filter(!is.na(p)) %>% arrange(p) %>% \n",
    "        transmute(phenotype, description, snp, OR = round(OR, \n",
    "            OR.digits), p = signif(p, p.digits), n_total, n_cases, \n",
    "            n_controls) %>% head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save our results\n",
    "\n",
    "We can save our results table and download it, as it does not contain participant-level information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write_csv(results_adj, file=\"adjusted_results.csv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.2"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}