meta_analysis/01_metal_meta_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perform the METAL meta-analysis\n",
    "\n",
    "In this notebook we perform the meta-analysis of aggregated *All of Us* and UK Biobank GWAS results.\n",
    "\n",
    "Note that this work is part of a larger project to [Demonstrate the Potential for Pooled Analysis of All of Us and UK Biobank Genomic Data](https://github.com/all-of-us/ukb-cross-analysis-demo-project). Specifically this is for the portion of the project that is the meta-analysis of **siloed** gwas results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "    <b>Cloud Environment</b>: This notebook was written for use on the <i>All of Us</i> Workbench.\n",
    "    <ul>\n",
    "        <li>Use \"Recommended Environment\" <kbd><b>General Analysis</b></kbd> which creates compute type <kbd><b>Standard VM</b></kbd> with reasonable defaults for CPU, RAM, and disk.</li>\n",
    "        <li>This notebook takes several minutes to run interactively. You can also it in the background via <kbd>run_notebook_in_the_background</kbd> for the sake of provenance and reproducibility.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "import os\n",
    "import pandas as pd\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = datetime.now()\n",
    "print(start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Install METAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Install METAL if it is not already installed.\n",
    "if [ ! -f ./generic-metal/metal ] ; then\n",
    "\n",
    "    curl -L -o metal.tar.gz \"http://csg.sph.umich.edu/abecasis/Metal/download/Linux-metal.tar.gz\"\n",
    "    tar -xf metal.tar.gz\n",
    "\n",
    "fi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Papermill parameters. See https://papermill.readthedocs.io/en/latest/usage-parameterize.html\n",
    "\n",
    "#---[ Inputs ]---\n",
    "AOU_GWAS_RESULTS = {\n",
    "    'HDL': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/aou/regenie/20230403/aou_alpha3_lipids_regenie_step2_HDL_norm.regenie',\n",
    "    'LDL': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/aou/regenie/20230403/aou_alpha3_lipids_regenie_step2_LDL_adjusted_norm.regenie',\n",
    "    'TC': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/aou/regenie/20230403/aou_alpha3_lipids_regenie_step2_TC_adjusted_norm.regenie',\n",
    "    'TG': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/aou/regenie/20230403/aou_alpha3_lipids_regenie_step2_TG_adjusted_norm.regenie'\n",
    "}\n",
    "\n",
    "UKB_GWAS_RESULTS = {\n",
    "    'HDL': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/ukb/regenie/20230403/ukb_200kwes_lipids_regenie_step2_HDL_mg_dl_norm.regenie',\n",
    "    'LDL': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/ukb/regenie/20230403/ukb_200kwes_lipids_regenie_step2_LDL_adj_mg_dl_norm.regenie',\n",
    "    'TC': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/ukb/regenie/20230403/ukb_200kwes_lipids_regenie_step2_TC_adj_mg_dl_norm.regenie',\n",
    "    'TG': 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/ukb/regenie/20230403/ukb_200kwes_lipids_regenie_step2_TG_log_mg_dl_norm.regenie'\n",
    "}\n",
    "\n",
    "METAL_RESULTS_BEFORE_INDEL_ID_FIXES = 'gs://fc-secure-e5c31994-13bb-4e3e-b7d7-db6effadc54f/data/metaanalysis/20230404/METAANALYSIS_LDL_1.tbl'\n",
    "\n",
    "#---[ Outputs ]---\n",
    "# Create a timestamp for a folder of results generated today.\n",
    "DATESTAMP = time.strftime('%Y%m%d')\n",
    "METAL_OUTPUTS = f'{os.getenv(\"WORKSPACE_BUCKET\")}/data/metaanalysis/{DATESTAMP}/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transfer inputs to local disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil -m cp {' '.join(AOU_GWAS_RESULTS.values())} ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil -m cp {' '.join(UKB_GWAS_RESULTS.values())} ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -lh *.regenie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fix ID column\n",
    "\n",
    "Fix UKB ids to use `chr_` format and also expand the ids for indels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head aou_alpha3_lipids_regenie_step2_HDL_norm.regenie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head ukb_200kwes_lipids_regenie_step2_HDL_mg_dl_norm.regenie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!grep ':I:' ukb_200kwes_lipids_regenie_step2_HDL_mg_dl_norm.regenie | head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!grep ':D:' ukb_200kwes_lipids_regenie_step2_HDL_mg_dl_norm.regenie | head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FIXED_ID_COL = 'FIXED_ID'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BASES = set(['A', 'C', 'G', 'T'])\n",
    "\n",
    "def expand_indel_ids(row):\n",
    "    try:\n",
    "        id_fields = row['ID'].split(':')\n",
    "        prefix = f'chr{id_fields[0]}_{id_fields[1]}'\n",
    "        \n",
    "        # ID holds bases for the ref allele.\n",
    "        if set(id_fields[2]).issubset(BASES):\n",
    "            assert set(id_fields[3]).issubset(BASES), f'{row[\"ID\"]}: {row[\"ALLELE0\"]} {row[\"ALLELE1\"]}'\n",
    "            # Just return ID in chr#_ format.\n",
    "            return f'{prefix}_{id_fields[2]}_{id_fields[3]}'\n",
    "        \n",
    "        # ID holds abbreviated notation for indels\n",
    "        change_len = round(float(id_fields[3])) # Some values were 3.01.\n",
    "        allele0_len = len(row['ALLELE0'])\n",
    "        allele1_len = len(row['ALLELE1'])\n",
    "        if 'I' == id_fields[2]:\n",
    "            if allele0_len < allele1_len:\n",
    "                ref = row['ALLELE0']\n",
    "                alt = row['ALLELE1']\n",
    "            elif allele0_len > allele1_len:\n",
    "                ref = row['ALLELE1']\n",
    "                alt = row['ALLELE0']\n",
    "            assert len(alt) - len(ref) == change_len, f'{row[\"ID\"]}: {row[\"ALLELE0\"]} {row[\"ALLELE1\"]}'\n",
    "        elif 'D' == id_fields[2]:\n",
    "            if allele0_len < allele1_len:\n",
    "                ref = row['ALLELE1']\n",
    "                alt = row['ALLELE0']\n",
    "            elif allele0_len > allele1_len:\n",
    "                ref = row['ALLELE0']\n",
    "                alt = row['ALLELE1']\n",
    "            assert len(ref) - len(alt) == change_len, f'{row[\"ID\"]}: {row[\"ALLELE0\"]} {row[\"ALLELE1\"]}'\n",
    "        else:\n",
    "            raise ValueError('unhandled id format')\n",
    "\n",
    "        # Return expanded form of ID in chr#_ format .\n",
    "        return f'{prefix}_{ref}_{alt}'\n",
    "    \n",
    "    except ValueError as e:\n",
    "        print(f'failed to fix id {row[\"ID\"]}: {row[\"ALLELE0\"]} {row[\"ALLELE1\"]}')\n",
    "        print(e)\n",
    "    except AssertionError as e:\n",
    "        print(f'failed to fix id {row[\"ID\"]}: {row[\"ALLELE0\"]} {row[\"ALLELE1\"]}')\n",
    "        print(e)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for gs_file in UKB_GWAS_RESULTS.values():\n",
    "    file = os.path.basename(gs_file)\n",
    "    df = pd.read_csv(file, sep=' ')\n",
    "    df[FIXED_ID_COL] = df.apply(expand_indel_ids, axis=1)\n",
    "    df.to_csv(file + '_' + FIXED_ID_COL, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head ukb_200kwes_lipids_regenie_step2_HDL_mg_dl_norm.regenie_{FIXED_ID_COL}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# METAL - Meta Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_metal(lipid):\n",
    "    aou_file = os.path.basename(AOU_GWAS_RESULTS[lipid])\n",
    "    # Use the file with the IDs in the 'chr' format.\n",
    "    ukb_file = os.path.basename(UKB_GWAS_RESULTS[lipid]) + '_' + FIXED_ID_COL\n",
    "    \n",
    "    \n",
    "    metal_parameters = f'''\n",
    "SCHEME STDERR\n",
    "AVERAGEFREQ ON\n",
    "MINMAXFREQ ON\n",
    "\n",
    "MARKER ID\n",
    "ALLELE ALLELE0 ALLELE1\n",
    "EFFECT BETA\n",
    "STDERR SE\n",
    "PVALUE Pvalue\n",
    "FREQ A1FREQ\n",
    "SEPARATOR WHITESPACE\n",
    "PROCESS {aou_file}\n",
    "\n",
    "MARKER {FIXED_ID_COL}\n",
    "ALLELE ALLELE0 ALLELE1\n",
    "EFFECT BETA\n",
    "STDERR SE\n",
    "PVALUE Pvalue\n",
    "FREQ A1FREQ\n",
    "SEPARATOR COMMA\n",
    "PROCESS {ukb_file}\n",
    "\n",
    "OUTFILE METAANALYSIS_{lipid}_ .tbl\n",
    "ANALYZE\n",
    "\n",
    "QUIT\n",
    "'''\n",
    "    print(f'Metal parameters:\\n{metal_parameters}')\n",
    "    \n",
    "    metal_parameters_filename = f'METAL_{lipid}.txt'\n",
    "    with open(metal_parameters_filename, 'w') as param_file:\n",
    "        param_file.write(metal_parameters)\n",
    "        \n",
    "    !./generic-metal/metal {metal_parameters_filename}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_metal('HDL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_metal('LDL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_metal('TC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_metal('TG')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls -lth METAANALYSIS_*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check that fixed ids are a subset of the unfixed ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_results = pd.read_table(METAL_RESULTS_BEFORE_INDEL_ID_FIXES)\n",
    "\n",
    "prior_results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "current_results = pd.read_table('METAANALYSIS_LDL_1.tbl')\n",
    "\n",
    "current_results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prior_ids = set(prior_results['MarkerName'].to_list())\n",
    "len(prior_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "current_ids = set(current_results['MarkerName'].to_list())\n",
    "len(current_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(prior_ids) > len(current_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newer_ids = current_ids - prior_ids\n",
    "len(newer_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newer_results = current_results[current_results['MarkerName'].isin(newer_ids)]\n",
    "\n",
    "newer_results.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newer_ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newer_results.groupby('Direction').count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'?+', '?-'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert set(newer_results['Direction'].to_list()) == set(['?+', '?-']), 'check that all \"new\" ids are from UKB only'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Store outputs in the workspace bucket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil -m cp METAANALYSIS_* {METAL_OUTPUTS}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gsutil ls -lh {METAL_OUTPUTS}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "end = datetime.now()\n",
    "print(end)\n",
    "print(end - start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Provenance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!date"
   ]
  }
 ],
 "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.12"
  },
  "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": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "353.266px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}