snps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b414298a",
   "metadata": {},
   "source": [
    "# SNP analysis\n",
    "After running GWAS on longevity, longevity with disease score covariances and all diseases, we run finemapping of the longevity with disease covariance.  \n",
    "Finemapping is performed via PolyFun.  \n",
    "Follow insuctions at: https://github.com/omerwe/polyfun  \n",
    "Finemapping results file was saved in output directory: output/finemap/polyfun_agg_longevity.txt.gz\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2a4f6ac7-d2df-4ff3-983b-847ebed6feb2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(rsid, allele1, allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(rsid, allele1, allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(rsid, allele1, allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(rsid, allele1, allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(rsid, allele1, allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(rsid, allele1, allele2)`\n"
     ]
    }
   ],
   "source": [
    "source(here::here(\"code/init.R\"))\n",
    "source(here::here(\"code/gwas.R\"))\n",
    "options(tgutil.cache=FALSE)\n",
    "pvals <- get_gwas_pvals() %cache_df% here(\"output/all_pvals.tsv\")\n",
    "snps <- pvals %>% filter(longevity_disease_covar_pval <=  log10(5e-8)) \n",
    "finemap_res <- tgutil::fread(here(\"output/finemap/polyfun_agg_longevity.txt.gz\")) %>% \n",
    "    separate(CREDIBLE_SET, c(\"chrom_locus\", \"start_locus\", \"end_locus\", \"CREDIBLE_SET\")) %>% \n",
    "    unite(chrom_locus:end_locus, col = \"locus\") %>% \n",
    "    as_tibble()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "125bef89",
   "metadata": {},
   "source": [
    "annotate finemapped snps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5c7f8a45",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, allele1, allele2)`\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A tibble: 11 × 2</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>chrom</th><th scope=col>n</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>chr1 </td><td>2</td></tr>\n",
       "\t<tr><td>chr14</td><td>1</td></tr>\n",
       "\t<tr><td>chr16</td><td>1</td></tr>\n",
       "\t<tr><td>chr18</td><td>1</td></tr>\n",
       "\t<tr><td>chr2 </td><td>1</td></tr>\n",
       "\t<tr><td>chr3 </td><td>1</td></tr>\n",
       "\t<tr><td>chr5 </td><td>2</td></tr>\n",
       "\t<tr><td>chr7 </td><td>1</td></tr>\n",
       "\t<tr><td>chr8 </td><td>1</td></tr>\n",
       "\t<tr><td>chr9 </td><td>1</td></tr>\n",
       "\t<tr><td>chrX </td><td>7</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 11 × 2\n",
       "\\begin{tabular}{ll}\n",
       " chrom & n\\\\\n",
       " <chr> & <int>\\\\\n",
       "\\hline\n",
       "\t chr1  & 2\\\\\n",
       "\t chr14 & 1\\\\\n",
       "\t chr16 & 1\\\\\n",
       "\t chr18 & 1\\\\\n",
       "\t chr2  & 1\\\\\n",
       "\t chr3  & 1\\\\\n",
       "\t chr5  & 2\\\\\n",
       "\t chr7  & 1\\\\\n",
       "\t chr8  & 1\\\\\n",
       "\t chr9  & 1\\\\\n",
       "\t chrX  & 7\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 11 × 2\n",
       "\n",
       "| chrom &lt;chr&gt; | n &lt;int&gt; |\n",
       "|---|---|\n",
       "| chr1  | 2 |\n",
       "| chr14 | 1 |\n",
       "| chr16 | 1 |\n",
       "| chr18 | 1 |\n",
       "| chr2  | 1 |\n",
       "| chr3  | 1 |\n",
       "| chr5  | 2 |\n",
       "| chr7  | 1 |\n",
       "| chr8  | 1 |\n",
       "| chr9  | 1 |\n",
       "| chrX  | 7 |\n",
       "\n"
      ],
      "text/plain": [
       "   chrom n\n",
       "1  chr1  2\n",
       "2  chr14 1\n",
       "3  chr16 1\n",
       "4  chr18 1\n",
       "5  chr2  1\n",
       "6  chr3  1\n",
       "7  chr5  2\n",
       "8  chr7  1\n",
       "9  chr8  1\n",
       "10 chr9  1\n",
       "11 chrX  7"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "19"
      ],
      "text/latex": [
       "19"
      ],
      "text/markdown": [
       "19"
      ],
      "text/plain": [
       "[1] 19"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "snps_annot <- snps %>%\n",
    "        left_join(finemap_res %>% mutate(chrom = paste0(\"chr\", CHR), start = BP, allele1 = A1, allele2 = A2)) %cache_df%\n",
    "        here(\"output/longevity_snps_finemapped.tsv\")\n",
    "snps_annot %>% \n",
    "    filter(is.na(PIP)) %>% \n",
    "    count(chrom)\n",
    "snps_annot %>% \n",
    "    filter(is.na(PIP)) %>% \n",
    "    nrow()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1006e1fd",
   "metadata": {},
   "source": [
    "Choose at least one locus from each ~1M region + high PIP loci:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "14245d98",
   "metadata": {},
   "outputs": [],
   "source": [
    "regs <- snps_annot %>%\n",
    "        gintervals.centers() %>%\n",
    "        gintervals.expand(5e5) %>%\n",
    "        gintervals.canonic()\n",
    "snps_top_annot <- snps_annot %>% \n",
    "    add_count(locus, name = \"n_locus\") %>%     \n",
    "    gintervals.neighbors1(regs) %>% \n",
    "    select(-dist) %>% \n",
    "    unite(\"region\", chrom1, start1, end1) %>% \n",
    "    add_count(region, name = \"n_region\") %>% \n",
    "    arrange(region, desc(PIP)) %>% \n",
    "    group_by(region) %>% \n",
    "    mutate(i = 1:n()) %>% \n",
    "    ungroup() %>% \n",
    "    filter(i == 1 | PIP >= 0.5) %>% \n",
    "    select(-i)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65bf6445",
   "metadata": {},
   "source": [
    "Compute number of significant SNPs 1MB around each top snp:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "315f3b8a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, allele1, allele2)`\n"
     ]
    }
   ],
   "source": [
    "snps_1M <- snps_top_annot %>%\n",
    "        mutate(orig_chrom = chrom, orig_start = start, orig_end = end) %>% \n",
    "        gintervals.centers() %>%\n",
    "        gintervals.expand(5e5) %>%\n",
    "        select(chrom, start, end, marker.ID, allele1, allele2, orig_chrom, orig_start, orig_end) %>%\n",
    "        gintervals.neighbors1(snps_annot %>% select(chrom, start, end), maxneighbors = 1e6) %>%\n",
    "        as_tibble() %>%\n",
    "        filter(dist == 0) %>%\n",
    "        count(chrom, start, end, marker.ID, allele1, allele2, orig_chrom, orig_start, orig_end, name = \"n_1M\") %>%\n",
    "        mutate(n_1M = n_1M - 1) %>%  # remove the SNP itself\n",
    "        select(-(chrom:end)) %>% \n",
    "        rename(chrom = orig_chrom, start = orig_start, end = orig_end)\n",
    "        \n",
    "snps_top_annot <- snps_top_annot %>% left_join(snps_1M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a0ba92e",
   "metadata": {},
   "source": [
    "Add results of COX regression (see GWAS_parents_survival notebook):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bba3a275",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, allele1, allele2,\n",
      "rsid)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, allele1, allele2,\n",
      "rsid)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, allele1, allele2,\n",
      "rsid)`\n"
     ]
    }
   ],
   "source": [
    "parents_mother <- fread(here(\"output/cox_parents_survival_mother_gwas.tsv\")) %>% \n",
    "    inner_join(snps_top_annot %>% select(chrom, start, end, rsid, marker.ID, allele1, allele2)) %>% \n",
    "    as_tibble()\n",
    "parents_father <- fread(here(\"output/cox_parents_survival_father_gwas.tsv\")) %>% \n",
    "    inner_join(snps_top_annot %>% select(chrom, start, end, rsid, marker.ID, allele1, allele2)) %>% \n",
    "    as_tibble()\n",
    "parents_both <- fread(here(\"output/cox_parents_survival_both_gwas.tsv\")) %>% \n",
    "    inner_join(snps_top_annot %>% select(chrom, start, end, rsid, marker.ID, allele1, allele2)) %>% \n",
    "    as_tibble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c09d650f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, rsid, allele1,\n",
      "allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, rsid, allele1,\n",
      "allele2)`\n",
      "\u001b[1m\u001b[22mJoining with `by = join_by(chrom, start, end, marker.ID, rsid, allele1,\n",
      "allele2)`\n"
     ]
    }
   ],
   "source": [
    "snps_top_annot <- snps_top_annot %>% \n",
    "    left_join(parents_mother %>% \n",
    "                select(chrom, start, end, rsid, marker.ID, allele1, allele2, cox_mother_surv_pval = pval, cox_mother_surv_z = z, cox_mother_surv_stat = Stat)             \n",
    "             ) %>% \n",
    "    left_join(parents_father %>% \n",
    "                select(chrom, start, end, rsid, marker.ID, allele1, allele2, cox_father_surv_pval = pval, cox_father_surv_z = z, cox_father_surv_stat = Stat)             \n",
    "             ) %>% \n",
    "    left_join(parents_both %>% \n",
    "                select(chrom, start, end, rsid, marker.ID, allele1, allele2, cox_both_surv_pval = pval, cox_both_surv_z = z, cox_both_surv_stat = Stat)             \n",
    "             ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cb9a4862",
   "metadata": {},
   "outputs": [],
   "source": [
    "fwrite(snps_top_annot, here(\"output/longevity_top_finemapped.csv\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fc96e2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "### H^2 SNP\n",
    "### With Diseases as covariates\n",
    "Create ldsc format sumstats:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65d9df49",
   "metadata": {},
   "source": [
    "adding std.err column from gwas of longevity with disease as covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "039dafef",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(marker.ID, allele1, allele2)`\n"
     ]
    }
   ],
   "source": [
    "gwas_longevity_disease_covar <- readr::read_rds(here(\"output/gwas_longevity_age_sex_disease_covar_extended.rds\"))\n",
    "pvals <- pvals %>% left_join(gwas_longevity_disease_covar %>% select(marker.ID, allele1, allele2, std.err))\n",
    "pvals %>%\n",
    "        filter(chrom != \"chrX\") %>% \n",
    "        mutate(N = 328542, CHR = gsub(\"chr\", \"\", chrom),  Z = longevity_disease_covar_beta / std.err) %>%\n",
    "        select(CHR, BP = start, A1 = allele1, A2 = allele2, N, Z, P = longevity_disease_covar_pval, SNP = rsid) %>%\n",
    "        fwrite(here(\"output/longevity_snps_ldsc.sumstats\"), sep = \" \", quote = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98662463",
   "metadata": {},
   "outputs": [],
   "source": [
    "at the terminal (polyfun conda):\n",
    "​\n",
    "./ldsc.py \\\n",
    "--out /home/nettam/projects/emr/ukbiobank/notebook/output/longevity_snps_ldsc.h2 \\\n",
    "--h2 /home/nettam/projects/emr/ukbiobank/notebook/output/longevity_snps_ldsc.sumstats \\\n",
    "--ref-ld-chr baselineLF2.2.UKB/baselineLF2.2.UKB. \\\n",
    "--w-ld-chr baselineLF2.2.UKB/weights.UKB. \\\n",
    "--not-M-5-50\n",
    "Total Observed scale h2: 0.1008 (0.0143)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02913c6c",
   "metadata": {},
   "source": [
    "### Without Diseases as covariates\n",
    "Create ldsc format sumstats: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "72898630",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[22mJoining with `by = join_by(marker.ID, allele1, allele2)`\n"
     ]
    }
   ],
   "source": [
    "gwas_longevity <- readr::read_rds(here(\"output/gwas_longevity_age_sex_covar_extended.rds\"))\n",
    "pvals <- pvals %>% left_join(gwas_longevity %>% select(marker.ID, allele1, allele2, std.err_longevity = std.err))\n",
    "pvals %>%\n",
    "        filter(chrom != \"chrX\") %>% \n",
    "        mutate(N = 328542, CHR = gsub(\"chr\", \"\", chrom),  Z = longevity_beta / std.err_longevity) %>%\n",
    "        select(CHR, BP = start, A1 = allele1, A2 = allele2, N, Z, P = longevity_pval, SNP = rsid) %>%\n",
    "        fwrite(here(\"output/longevity_snps_no_covar_ldsc.sumstats\"), sep = \" \", quote = FALSE)\n",
    "     \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "656fd052",
   "metadata": {},
   "outputs": [],
   "source": [
    "at the terminal (polyfun conda):\n",
    "\n",
    "./ldsc.py \\\n",
    "--out /home/nettam/projects/emr/ukbiobank/notebook/output/longevity_snps_no_covar_ldsc.h2 \\\n",
    "--h2 /home/nettam/projects/emr/ukbiobank/notebook/output/longevity_snps_no_covar_ldsc.sumstats \\\n",
    "--ref-ld-chr baselineLF2.2.UKB/baselineLF2.2.UKB. \\\n",
    "--w-ld-chr baselineLF2.2.UKB/weights.UKB. \\\n",
    "--not-M-5-50\n",
    "Total Observed scale h2: 0.1641 (0.0139)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 4.2",
   "language": "R",
   "name": "ir42"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}