workflow/gwas_enrichment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "# Enrichment analysis workflows\n",
    "\n",
    "This workflow document has several pipelines, all written in SoS Workflow language. The `torus` pipeline is based on some `snakemake` pipelines originally written by Jean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <table class=\"revision_table\">\n",
       "        <tr>\n",
       "        <th>Revision</th>\n",
       "        <th>Author</th>\n",
       "        <th>Date</th>\n",
       "        <th>Message</th>\n",
       "        <tr>\n",
       "        <tr><td><a target=\"_blank\" href=\"git@github.com:gaow/fine-mapping/blob/4e3cda48d32cd317d5c62ae037051e5f4f8c2659/workflow/gwas_enrichment.ipynb\"><span class=\"revision_id\">4e3cda4<span></a></td>\n",
       "<td>minqiao</td>\n",
       "<td>2019-07-09</td>\n",
       "<td>update pipeline</td></tr><tr><td><a target=\"_blank\" href=\"git@github.com:gaow/fine-mapping/blob/8200a370c7292ff6830a9b60e6b41341846b09d1/workflow/gwas_enrichment.ipynb\"><span class=\"revision_id\">8200a37<span></a></td>\n",
       "<td>Gao Wang</td>\n",
       "<td>2019-07-07</td>\n",
       "<td>Add enrichment workflow</td></tr></table>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%revisions -s -n 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: sos run gwas_enrichment.ipynb\n",
      "               [workflow_name | -t targets] [options] [workflow_options]\n",
      "  workflow_name:        Single or combined workflows defined in this script\n",
      "  targets:              One or more targets to generate\n",
      "  options:              Single-hyphen sos parameters (see \"sos run -h\" for details)\n",
      "  workflow_options:     Double-hyphen workflow-specific parameters\n",
      "\n",
      "Workflows:\n",
      "  zscore2bed\n",
      "  get_variants\n",
      "  merge_annotations\n",
      "  range2var_annotation\n",
      "  enrichment\n",
      "  deepsea_apply\n",
      "  gregor\n",
      "\n",
      "Global Workflow Options:\n",
      "  --deepsea-model  path(f'{cwd}/deepsea/deepsea_ATAC053018.hdf5')\n",
      "\n",
      "                        Deepsea model for use with `deepsee_apply`\n",
      "  --annotation-dir  f'{cwd}/bed_annotations'\n",
      "\n",
      "                        Path to directory of annotation files\n",
      "  --z-score . (as path)\n",
      "                        Path to z-score file\n",
      "  --single-annot . (as path)\n",
      "                        Path to list of single annotations to use\n",
      "  --multi-annot  paths() #parameter: multi_annot = paths([\"data/all_annotations.txt\", \"data/general_annotations.txt\", \"data/reduced_annotations.txt\"])\n",
      "\n",
      "                        Path to lists of multi-annotations to use\n",
      "\n",
      "Sections\n",
      "  zscore2bed_1:         Auxiliary step to get variant in bed format based on\n",
      "                        variant ID in z-score file\n",
      "    Workflow Options:\n",
      "      --in-file . (as path)\n",
      "  zscore2bed_2:\n",
      "  get_variants:\n",
      "  merge_annotations:    Auxiliary step to merge multiple annotations\n",
      "    Workflow Options:\n",
      "      --out-file . (as path)\n",
      "      --data-files paths()\n",
      "\n",
      "  range2var_annotation_1: Get variants in data that falls in target region\n",
      "  range2var_annotation_2: Make binary annotation file\n",
      "    Workflow Options:\n",
      "      --discrete 1 (as int)\n",
      "  range2var_annotation_3: Obtain multiple annotations per SNP for enrichment\n",
      "                        analysis\n",
      "  enrichment_1:         Run torus with annotation + z-score file\n",
      "    Workflow Options:\n",
      "      --gmap  path(f\"{cwd}/gene_map/gmap.gz\") \n",
      "\n",
      "  enrichment_2:         Consolidate all torus result into one table\n",
      "  deepsea_apply_1:\n",
      "    Workflow Options:\n",
      "      --variants . (as path)\n",
      "  deepsea_apply_2:\n",
      "  gregor_1:\n",
      "    Workflow Options:\n",
      "      --gregor-input . (as path)\n",
      "  gregor_2:\n",
      "    Workflow Options:\n",
      "      --gregor-db /home/min/Documents/hg19/GREGOR_DB (as path)\n",
      "      --pop EUR\n",
      "  gregor_3:\n",
      "    Workflow Options:\n",
      "      --gregor-path /home/min/Documents/GREGOR (as path)\n",
      "  gregor_4:\n"
     ]
    }
   ],
   "source": [
    "!sos run gwas_enrichment.ipynb -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Copy this GIT depository\n",
    "```\n",
    "git clone https://github.com/gaow/fine-mapping.git\n",
    "```\n",
    "\n",
    "## To run enrichment analysis (overview)\n",
    "\n",
    "### Step 1 (start from `bed` format annotation)\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb range2var_annotation --z-score ... --single-annot ... --multi-annot ...\n",
    "```\n",
    "### Step 2 (start from binary annotation for SNPs)\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb enrichment --z-score ... --single-annot ... --multi-annot ...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Software intallation\n",
    "\n",
    "### `SoS`\n",
    "    \n",
    "SoS introduction [webpage](https://vatlab.github.io/sos-docs/)\n",
    "\n",
    "Installation:\n",
    "\n",
    "```\n",
    "pip install sos\n",
    "```\n",
    "\n",
    "SoS basic usage [webpage](https://vatlab.github.io/sos-docs/running.html#content)\n",
    "\n",
    "### `torus`\n",
    "\n",
    "`torus` introduction and example [webpage](https://github.com/xqwen/torus/)\n",
    "    \n",
    "Installation:\n",
    "\n",
    "```\n",
    "git clone https://github.com/xqwen/torus.git\n",
    "cd torus/src/\n",
    "make\n",
    "chmod +x torus\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Reference data preparation\n",
    "### hg19.fa\n",
    "- only for `deepsea` and `gregor` steps; if you only use enrichment analysis, you do not need to download it.\n",
    "- Genome Reference: hg19 (GRCh37) or hg38 (GRCh38). \n",
    "- Download [link](http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Data specification\n",
    "### Annotation specification\n",
    "#### `bed` format annotation input\n",
    "**If you use this type of annotation file, you need to run step 1 and 2.**\n",
    "\n",
    "Annotation files are in `bed` format (for example: Promoter_UCSC.bed).\n",
    "\n",
    "        chr1\t9873\t16361\n",
    "        chr1\t32610\t36610\n",
    "        chr1\t67090\t71090\n",
    "        ...\n",
    "\n",
    "#### Binary annotations of interest\n",
    "**If you use this type of annotation file, you only need to run step 2.**\n",
    "\n",
    "0/1 binary indicated. For example, a few lines in `Coding_UCSC_annotation.torus.gz`.\n",
    "\n",
    "    SNP     Coding_UCSC_d\n",
    "    9:140192349:G:C    0\n",
    "    9:140192576:C:T    0\n",
    "    9:140194020:G:A    0\n",
    "    9:140194735:G:T    1\n",
    "    9:140194793:A:G    1\n",
    "    9:140195019:C:A    1\n",
    "\n",
    "\n",
    "There are many annotations one can use. For this workflow you should prepare two lists:\n",
    "\n",
    "- For the ones that you'd like to use for one variable `torus`, put their file name prefix in `data/single_annotations.txt`\n",
    "- For the ones that you'd like to use for multi-variable `torus`, put their file name prefix in `data/multi_annotations.txt`\n",
    "\n",
    "If you want to remove one annotation just use `#` to comment it out in the list files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div class=\"sos_hint\">> ../data/general_annotations.txt (100 B):</div>"
      ],
      "text/plain": [
       "\n",
       "> ../data/general_annotations.txt (100 B):"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div class=\"sos_hint\">6 lines</div>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coding_UCSC\n",
      "Promoter_UCSC\n",
      "Repressed_Hoffman\n",
      "Conserved_LindbladToh\n",
      "DHS_peaks_Trynka\n",
      "#FetalDHS_Trynka"
     ]
    }
   ],
   "source": [
    "%preview ../data/general_annotations.txt -l 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "### GWAS z-scores specification\n",
    "\n",
    "#### GWAS data of interest\n",
    "\n",
    "Example GWAS data (European)\n",
    "\n",
    "[SCZ Sweden data](https://www.med.unc.edu/pgc/results-and-downloads/data-use-agreement-forms/SwedenSCZ_data_download_agreement)\n",
    "\n",
    "[PTSD European Ancestry](https://www.med.unc.edu/pgc/results-and-downloads/data-use-agreement-forms/PTSD%20EA_data_download_agreement)\n",
    "\n",
    "#### LD block files\n",
    "\n",
    "The last column indicates LD chunk. LD chunk blocks depend on chromosome and position. There are 1703 LD chunks in European human genome.\n",
    "\n",
    "```\n",
    "            chr1    10583       1892607     1\n",
    "            chr1    1892607     3582736     2\n",
    "            ...\n",
    "            chr22   49824534    51243298    1703\n",
    "```\n",
    "\n",
    "#### GWAS input format\n",
    "\n",
    "Format `chr:pos:alt:ref   ld_label   z-score`. 3 columns in total: `chr:pos:alt:ref` (chromosome, position, alternative allele and reference allele), `ld_label` (LD chunk label) and `z-score` (summary statistics).\n",
    "\n",
    "Example data:\n",
    "```\n",
    "1:10583:A:G        1      0.116319\n",
    "...\n",
    "22:51228910:A:G    1703   -0.866894\n",
    "```\n",
    "\n",
    "All GWAS summary statistics have to be converted to this format in order to use.\n",
    "\n",
    "### SNP map\n",
    "GWAS SNP map. Identical number of SNPs with GWAS.\n",
    "\n",
    "Format `chrom.pos  chr  pos`. 3 columns in total: `chrom.pos` (example: `chr1.729679`), `chr` and `pos`.\n",
    "```\n",
    "chr1.729679\t1\t729679\n",
    "chr1.731718\t1\t731718\n",
    "chr1.734349\t1\t734349\n",
    "...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Minimal working example file list\n",
    "### GWAS\n",
    "GWAS example [link](https://www.dropbox.com/s/y2q29a71bb2wmr1/SCZGWAS.zscore.gz?dl=0)\n",
    "```\n",
    "1:729679:C:G\t1\t-1.2583\n",
    "1:731718:T:C\t1\t0.9745\n",
    "1:734349:T:C\t1\t1.0247\n",
    "...\n",
    "```\n",
    "\n",
    "### SNP map\n",
    "SNP map example [link](https://www.dropbox.com/s/oi4yr0au3oy3hi2/smap.gz?dl=0)\n",
    "```\n",
    "chr1.729679\t1\t729679\n",
    "chr1.731718\t1\t731718\n",
    "chr1.734349\t1\t734349\n",
    "...\n",
    "```\n",
    "\n",
    "### Functional annotation\n",
    "Example: `Promoter_UCSC.bed` download [link](https://www.dropbox.com/s/huo5k78wptxmrpe/Promoter_UCSC.bed?dl=0). 22,436 lines in total. First few lines are as follows:\n",
    "\n",
    "        chr1\t9873\t16361\n",
    "        chr1\t32610\t36610\n",
    "        chr1\t67090\t71090\n",
    "        ...\n",
    "\n",
    "### LD blocks\n",
    "Download [link](https://www.dropbox.com/s/up78cnpkpoodark/ld_chunk.bed?dl=0). Used to identify LD chunk for each GWAS SNP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Workflow\n",
    "Test on my desktop\n",
    "```\n",
    "cd ~/GIT/fine-mapping\n",
    "sos run workflow/gwas_enrichment.ipynb range2var_annotation --cwd ~/Documents/GWAS_ATAC --annotation_dir ~/Documents/GWAS_ATAC/bed_annotations \\\n",
    "--z-score ~/Documents/GWAS_ATAC/SCZGWAS_zscore/SCZGWAS.zscore.gz --single-annot ~/GIT/fine-mapping/data/general_annotations.txt -v3\n",
    "sos run workflow/gwas_enrichment.ipynb enrichment --cwd ~/Documents/GWAS_ATAC --annotation_dir ~/Documents/GWAS_ATAC/bed_annotations \\\n",
    "--z-score ~/Documents/GWAS_ATAC/SCZGWAS_zscore/SCZGWAS.zscore.gz --single-annot ~/GIT/fine-mapping/data/general_annotations.txt -v3\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[global]\n",
    "# working directory\n",
    "parameter: cwd = path() # 1\n",
    "# hg19 path\n",
    "parameter: hg19 = path()\n",
    "# Deepsea model for use with `deepsee_apply`\n",
    "parameter: deepsea_model = path()\n",
    "# Path to directory of annotation files\n",
    "parameter: annotation_dir = path() # 2\n",
    "# Path to z-score file\n",
    "parameter: z_score = path() # 3\n",
    "# Path to list of single annotations to use\n",
    "parameter: single_annot = path() # 4\n",
    "# Path to lists of multi-annotations to use\n",
    "parameter: multi_annot = paths() # 5\n",
    "try:\n",
    "    single_anno = [f\"{annotation_dir}/{x.split()[0]}.bed\" for x in open(single_annot).readlines() if not x.startswith('#')]\n",
    "    multi_anno = dict([(f'{y:bn}', [f\"{annotation_dir}/{x.split()[0]}.bed\" for x in open(y).readlines() if not x.startswith('#')]) for y in multi_annot])\n",
    "except (FileNotFoundError, IsADirectoryError):\n",
    "    single_anno = []\n",
    "    multi_anno = dict()\n",
    "out_dir = f'{cwd}/{z_score:bn}'.replace('.', '_')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Utility steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "### Convert variants from z-score file to bed format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Auxiliary step to get variant in bed format based on variant ID in z-score file\n",
    "[zscore2bed_1]\n",
    "depends: R_library('readr'), R_library('stringr'), R_library('dplyr>=0.7.7')\n",
    "parameter: in_file = path()\n",
    "input: in_file\n",
    "output: f'{_input:n}.bed.unsorted'\n",
    "R: expand = \"${ }\", container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout'\n",
    "    library(readr)\n",
    "    library(stringr)\n",
    "    library(dplyr)\n",
    "    var_file <- ${_input:r}\n",
    "    out_file <- ${_output:r}\n",
    "\n",
    "    variants <- read_tsv(var_file, col_names=FALSE)\n",
    "    colnames(variants) = c('variant', 'ld', 'zscore')\n",
    "    var_info <- str_split(variants$variant, \":\")\n",
    "    variants <- mutate(variants, chr = paste0(\"chr\", sapply(var_info, function(x){x[1]})), \n",
    "                                 pos = sapply(var_info, function(x){x[2]})) %>%\n",
    "                mutate(start = as.numeric(pos), stop=as.numeric(pos) + 1) %>%\n",
    "                select(chr, start, stop, variant)\n",
    "    options(scipen=1000) # So that positions are always fully written out)\n",
    "    write.table(variants, file=out_file, quote=FALSE, col.names=FALSE, row.names=FALSE, sep=\"\\t\")\n",
    "\n",
    "[zscore2bed_2]\n",
    "depends: executable('sort-bed')\n",
    "output: f'{_input:n}'\n",
    "bash: expand = True, container = 'gaow/atac-gwas', workdir = cwd\n",
    "     sort-bed {_input} > {_output}\n",
    "\n",
    "[get_variants: provides = '{data}.bed']\n",
    "output: f'{data}.bed'\n",
    "sos_run('zscore2bed', in_file = f'{_output:n}.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "### Merge annotations for joint torus analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Auxiliary step to merge multiple annotations\n",
    "[merge_annotations]\n",
    "depends: R_library('readr'), R_library('dplyr>=0.7.7'), R_library('stringr')\n",
    "parameter: out_file=path()\n",
    "parameter: data_files=paths()\n",
    "input: data_files\n",
    "output: out_file\n",
    "R: expand = '${ }', container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout'\n",
    "\n",
    "    library(readr)\n",
    "    library(dplyr)\n",
    "    library(stringr)\n",
    "\n",
    "    out_name <- ${_output:r}\n",
    "    annots <- c(${_input:r,})\n",
    "    var_annot <- read_tsv(annots[1], col_names=T)\n",
    "    for(i in 2:length(annots)){\n",
    "        var_annot2 <- read_tsv(annots[i], col_names=T)\n",
    "        stopifnot(all(var_annot$SNP == var_annot2$SNP))\n",
    "        var_annot <- cbind(var_annot, var_annot2[,2])\n",
    "    }\n",
    "    write.table(var_annot, file=gzfile(out_name),\n",
    "                row.names=FALSE, quote=FALSE, sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Prepare `torus` format binary annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Get variants in data that falls in target region\n",
    "[range2var_annotation_1]\n",
    "depends: f'{z_score:n}.bed', executable('bedops')\n",
    "input: set(paths(single_anno + list(multi_anno.values()))), group_by = 1\n",
    "output: f'{out_dir}/{_input:bn}.{z_score:bn}.bed'\n",
    "bash: expand = True, container = 'gaow/atac-gwas', workdir = cwd\n",
    "    bedops -e {z_score:n}.bed {_input} > {_output}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Make binary annotation file\n",
    "[range2var_annotation_2]\n",
    "parameter: discrete = 1\n",
    "depends: z_score, R_library('readr'), R_library('dplyr>=0.7.7'), R_library('stringr')\n",
    "input: \n",
    "output: f'{_input:n}.gz'\n",
    "R: expand = \"${ }\", container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout'\n",
    "    library(readr)\n",
    "    library(dplyr)\n",
    "    library(stringr)\n",
    "\n",
    "    variant_tsv <- ${z_score:r}\n",
    "    annotation_var_bed <- ${_input:r}\n",
    "    annot_name <- ${_input:bnr} %>% str_replace(paste0(\".\",${z_score:bnr}), \"\")\n",
    "    out_name <- ${_output:r}\n",
    "\n",
    "    vars <- read_tsv(variant_tsv, col_names=FALSE)[,1]\n",
    "    annot_vars = read_tsv(annotation_var_bed, col_names=FALSE)\n",
    "    names(vars) <- \"SNP\"\n",
    "    vars <- vars %>%\n",
    "            mutate(annot_d = case_when(SNP %in% annot_vars$X4 ~ 1,\n",
    "                                                        TRUE ~ 0))\n",
    "    names(vars)[2] <- paste0(annot_name, '${\"_d\" if discrete else \"_c\"}')\n",
    "    write.table(vars, file=gzfile(out_name),\n",
    "                col.names=TRUE, row.names=FALSE, sep=\"\\t\", quote=FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Obtain multiple annotations per SNP for enrichment analysis\n",
    "[range2var_annotation_3]\n",
    "for k, value in multi_anno.items():\n",
    "    sos_run('merge_annotations', out_file = f'{out_dir}/{k}.{z_score:bn}.gz', data_files = [f'{out_dir}/{v:bn}.{z_score:bn}.gz' for v in paths(value)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Enrichment analysis via `torus`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Run torus with annotation + z-score file\n",
    "[enrichment_1 (run torus)]\n",
    "depends: z_score\n",
    "parameter: gmap = path(f\"{cwd}/gene_map/gmap.gz\") \n",
    "input_files = [f'{out_dir}/{value:bn}.{z_score:bn}.gz' for value in paths(single_anno)] + [f'{out_dir}/{k}.{z_score:bn}.gz' for k in multi_anno]\n",
    "fail_if(len(input_files) == 0, msg = \"No annotations to use! Please use ``--single-annot`` and ``--multi-annot`` to pass annotation lists\")\n",
    "input: input_files, group_by = 1\n",
    "output: f'{_input:n}.torus'\n",
    "\n",
    "bash: container = 'gaow/atac-gwas', workdir = cwd, expand = True, stdout = f'{_output:n}.stdout', stderr = f'{_output:n}.stderr'\n",
    "    rm -rf {_output:n}_prior \n",
    "    torus -d {z_score} -smap {_input:d}/smap.gz -gmap {gmap} -annot {_input} -est --load_zval -dump_prior {_output:n}_prior > {_output}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "# Consolidate all torus result into one table\n",
    "[enrichment_2 (make output table)]\n",
    "depends: R_library('readr'), R_library('dplyr>=0.7.7'), R_library('stringr')\n",
    "input: group_by = 'all'\n",
    "output: f'{out_dir}/{z_score:bn}.torus.merged.gz'\n",
    "R: expand = '${ }', container = 'gaow/atac-gwas', workdir = cwd, stdout = f'{_output:n}.stdout'\n",
    "    library(dplyr)\n",
    "    library(readr)\n",
    "    library(stringr)\n",
    "    output_file <- ${_output:r}\n",
    "    names <- c(${_input:bnr,})\n",
    "    file_names <- c(${_input:r,})\n",
    "    tab <- data.frame(\"name\"=names)\n",
    "\n",
    "    data <- lapply(file_names, function(f){\n",
    "                    lns <- read_lines(f)\n",
    "                    x <- unlist(str_split(lns[3], \" \"))\n",
    "                    x <- x[str_length(x) > 0]\n",
    "                    x})\n",
    "\n",
    "    tab[,\"log odds ratio\"] <- sapply(data, function(x){x[2]})\n",
    "    tab[,\"95% ci\"] <- sapply(data, function(x){paste0(\"(\", x[3], \", \", x[4], \")\")})\n",
    "    write.table(tab, file=gzfile(output_file), sep=\"\\t\", col.names=TRUE, row.names=FALSE, quote=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Apply pre-trained `deepsea` model to variants\n",
    "\n",
    "Credits to Yanyu Liang. Required inputs are:\n",
    "\n",
    "- Path to hg19 reference genome\n",
    "- Path to model HDF5 file\n",
    "- Path to list of variants, in the format of:\n",
    "\n",
    "```\n",
    "chr1\t68090\t68091\tG\tT\n",
    "chr1\t68090\t68091\tG\tA\n",
    "chr1\t68090\t68091\tG\tC\n",
    "chr1\t68091\t68092\tC\tG\n",
    "chr1\t68091\t68092\tC\tT\n",
    "chr1\t68091\t68092\tC\tA\n",
    "chr1\t68092\t68093\tT\tA\n",
    "chr1\t68092\t68093\tT\tC\n",
    "chr1\t68092\t68093\tT\tG\n",
    "chr1\t68093\t68094\tA\tC\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "To run this workflow:\n",
    "    \n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants /path/to/variant/list.file\n",
    "```\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/scz.swe.var.txt\n",
    "```\n",
    "It requires much memory. So we test for first 10K line of variant and it worked.\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/test.var.txt\n",
    "sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/test100K.var.txt\n",
    "```\n",
    "Break 9898079 snps into two halves and run deepsea separately.\n",
    "```\n",
    "less scz.swe.var.txt | head -4949039 | gzip > half1.var.txt\n",
    "```\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/half1.var.txt\n",
    "sos run workflow/gwas_enrichment.ipynb deepsea_apply --variants ~/Documents/GWAS_ATAC/GWAS_data/SCZSweden/half2.var.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[deepsea_apply_1 (prepare config file)]\n",
    "parameter: variants = path()\n",
    "input: variants\n",
    "output: f'{_input:n}.config'\n",
    "report: expand = True, output = _output\n",
    "    var_list: {_input:r}\n",
    "    pred_model:\n",
    "      path: '{deepsea_model}'\n",
    "      label:\n",
    "        CN: 1\n",
    "        DN: 2\n",
    "        GA: 3\n",
    "        ips: 4\n",
    "        NSC: 5\n",
    "    reference_genome: '{hg19}'\n",
    "    script_dir: '/opt/deepann'\n",
    "    out_dir: {_output:dr}\n",
    "    out_prefix: '{_output:bn}_deepsea'\n",
    "\n",
    "    ##### some default setup #####\n",
    "    #### usually don't change ####\n",
    "    check_allele: False\n",
    "    design: '499-1-500'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[deepsea_apply_2 (apply deepsea weights)]\n",
    "depends: executable('snakemake')\n",
    "output: f'{_input:n}.deepsea.list'\n",
    "bash: expand = True, container = 'gaow/deepann', volumes = [f'{hg19:d}:{hg19:d}', f'{deepsea_model:d}:{deepsea_model:d}', f'{os.path.expanduser(\"~\")}:/home/$USER'], extra_args = \"-e HOME=/home/$USER -e USER=$USER\"\n",
    "    snakemake --snakefile /opt/deepann/Snakefile --configfile {_input} && ls {_input:n}_deepsea/output__*.gz > {_output}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Enrichment analysis via GREGOR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "To properly perform enrichment analysis we want to match the control SNPs with the SNPs of interest -- that is, SNPs inside CS -- in terms of LD, distance to nearest gene and MAF. The [GREGOR](http://csg.sph.umich.edu/GREGOR/index.php/site/download) software can generate list of matched SNPs. I will use SNPs inside CS as input and expect a list of output SNPs matching these inputs.\n",
    "\n",
    "GREGOR is release under University of Michigan license so I'll not make it into a docker image. So path to GREGOR directory is required. Also we need reference files, prepared by:\n",
    "\n",
    "```\n",
    "cat \\\n",
    "    GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz.part.00 \\\n",
    "    GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz.part.01 \\\n",
    "    > GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz\n",
    "tar zxvf GREGOR.AFR.ref.r2.greater.than.0.7.tar.gz\n",
    "```\n",
    "\n",
    "MD5SUM check:\n",
    "\n",
    "```\n",
    "AFR.part.0\t( MD5: 9926904128dd58d6bf1ad4f1e90638af )\n",
    "AFR.part.1\t( MD5: c1d30aff89a584bfa8c1fa1bdc197f21 )\n",
    "```\n",
    "\n",
    "Same for EUR data-set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "#### To run GREGOR\n",
    "`PGC3`\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb gregor --gregor-input ~/Documents/GWAS_ATAC/scz2_zscore/gregor.scz2.txt.gz --single-annot data/gregor.txt\n",
    "```\n",
    "`\n",
    "PGC2\n",
    "`\n",
    "```\n",
    "sos run workflow/gwas_enrichment.ipynb gregor --gregor-input ~/Documents/GWAS_ATAC/pgc2_zscore/gregor.pgc2.txt.gz --single-annot data/gregor.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[gregor_1 (make SNP index)]\n",
    "depends: executable('zcat')\n",
    "parameter: gregor_input = path()\n",
    "input: gregor_input\n",
    "output: f'{_input:nn}.rsid.txt', f'{_input:nn}.annotations.list'\n",
    "bash: expand = '${ }'\n",
    "    zcat ${_input} | cut -f 2,3 -d \"_\" | sed 's/_/:/g' > ${_output[0]}\n",
    "\n",
    "with open(_output[1], 'w') as f:\n",
    "    f.write('\\n'.join(single_anno))\n",
    "\n",
    "[gregor_2 (make configuration file)]\n",
    "depends: executable('sed')\n",
    "parameter: gregor_db = path('~/Documents/hg19/GREGOR_DB')\n",
    "parameter: pop = 'EUR'\n",
    "output: f'{_input[0]:nn}.gregor.conf'\n",
    "report: output = f'{_output}', expand = True\n",
    "    ##############################################################################\n",
    "    # CHIPSEQ ENRICHMENT CONFIGURATION FILE\n",
    "    # This configuration file contains run-time configuration of\n",
    "    # CHIP_SEQ ENRICHMENT\n",
    "    ###############################################################################\n",
    "    ## KEY ELEMENTS TO CONFIGURE : NEED TO MODIFY\n",
    "    ###############################################################################\n",
    "    INDEX_SNP_FILE = {_input[0]}\n",
    "    BED_FILE_INDEX = {_input[1]} \n",
    "    REF_DIR = {gregor_db}\n",
    "    R2THRESHOLD = 0.7 ## must be greater than 0.7\n",
    "    LDWINDOWSIZE = 10000 ## must be less than 1MB; these two values define LD buddies\n",
    "    OUT_DIR = {_output:nn}_gregor_output\n",
    "    MIN_NEIGHBOR_NUM = 10 ## define the size of neighborhood\n",
    "    BEDFILE_IS_SORTED = true  ## false, if the bed files are not sorted\n",
    "    POPULATION = {pop}  ## define the population, you can specify EUR, AFR, AMR or ASN\n",
    "    TOPNBEDFILES = 2 \n",
    "    JOBNUMBER = 10\n",
    "    ###############################################################################\n",
    "    #BATCHTYPE = mosix ##  submit jobs on MOSIX\n",
    "    #BATCHOPTS = -E/tmp -i -m2000 -j10,11,12,13,14,15,16,17,18,19,120,122,123,124,125 sh -c\n",
    "    ###############################################################################\n",
    "    #BATCHTYPE = slurm   ##  submit jobs on SLURM\n",
    "    #BATCHOPTS = --partition=broadwl --account=pi-mstephens --time=0:30:0\n",
    "    ###############################################################################\n",
    "    BATCHTYPE = local ##  run jobs on local machine\n",
    "\n",
    "bash: expand = True\n",
    "    sed -i '/^$/d' {_output}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "GREGOR is written in `perl`. Some libraries are required before one can run GREGOR:\n",
    "\n",
    "```\n",
    "sudo apt-get install libdbi-perl libswitch-perl libdbd-sqlite3-perl\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[gregor_3 (run gregor)]\n",
    "depends: executable('perl')\n",
    "parameter: gregor_path = path('~/Documents/GREGOR')\n",
    "output: f'{_input:nn}_gregor_output/StatisticSummaryFile.txt'\n",
    "bash: expand = True\n",
    "    perl {gregor_path}/script/GREGOR.pl --conf {_input} && touch {_output}\n",
    "\n",
    "[gregor_4 (format output)]\n",
    "depends: executable('sed')\n",
    "output: f'{_input:n}.csv'\n",
    "bash: expand = True\n",
    "    sed 's/\\t/,/g' {_input} > {_output}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SoS",
   "language": "sos",
   "name": "sos"
  },
  "language_info": {
   "codemirror_mode": "sos",
   "file_extension": ".sos",
   "mimetype": "text/x-sos",
   "name": "sos",
   "nbconvert_exporter": "sos_notebook.converter.SoS_Exporter",
   "pygments_lexer": "sos"
  },
  "sos": {
   "default_kernel": "SoS",
   "kernels": [
    [
     "Python3",
     "python3",
     "Python3",
     "#FFD91A"
    ],
    [
     "SoS",
     "sos",
     "",
     ""
    ]
   ],
   "panel": {
    "displayed": false,
    "height": 0,
    "style": "side"
   },
   "version": "0.19.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}