{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9e70a797-4c24-4c7c-bae5-f2b94b16dbd3",
   "metadata": {
    "kernel": "SoS",
    "tags": []
   },
   "source": [
    "# Other colocalization methods\n",
    "\n",
    "## Goal\n",
    "\n",
    "This notebook takes the input and run Hyprcoloc, MOLOC or COLOC(V5) (SuSiE-Coloc) , the result can be summarized to get power and FDR.\n",
    "\n",
    "Because different methods have different output, we summarized four elements to calculate FDR and power - coloc_trait, coloc_set true_trait and true_variant.\n",
    "\n",
    "## Input \n",
    "\n",
    "Individual level data X and Y or summary statistics z and LD from other notebooks. put in the parameter: simufile.\n",
    "\n",
    "## Output\n",
    "\n",
    "Hyprcoloc / MOLOC / COLOC (V5) original result, along with some summarized elements."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b60160a4-a472-44c2-80d8-47d430a2337c",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "Example output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4f48e1ca-bd4a-4cb3-8f74-5fcf202b0610",
   "metadata": {
    "kernel": "R",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 3 × 7</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>iteration</th><th scope=col>traits</th><th scope=col>posterior_prob</th><th scope=col>regional_prob</th><th scope=col>candidate_snp</th><th scope=col>posterior_explained_by_snp</th><th scope=col>dropped_trait</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>1</td><td>5, 9</td><td>0.8754</td><td>1</td><td>snp917 </td><td>0.3881</td><td>NA</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>2</td><td>3, 8</td><td>0.7915</td><td>1</td><td>snp2568</td><td>0.0418</td><td>NA</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>3</td><td>2, 6</td><td>0.9426</td><td>1</td><td>snp2904</td><td>0.1141</td><td>NA</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 3 × 7\n",
       "\\begin{tabular}{r|lllllll}\n",
       "  & iteration & traits & posterior\\_prob & regional\\_prob & candidate\\_snp & posterior\\_explained\\_by\\_snp & dropped\\_trait\\\\\n",
       "  & <dbl> & <chr> & <dbl> & <dbl> & <chr> & <dbl> & <chr>\\\\\n",
       "\\hline\n",
       "\t1 & 1 & 5, 9 & 0.8754 & 1 & snp917  & 0.3881 & NA\\\\\n",
       "\t2 & 2 & 3, 8 & 0.7915 & 1 & snp2568 & 0.0418 & NA\\\\\n",
       "\t3 & 3 & 2, 6 & 0.9426 & 1 & snp2904 & 0.1141 & NA\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 3 × 7\n",
       "\n",
       "| <!--/--> | iteration &lt;dbl&gt; | traits &lt;chr&gt; | posterior_prob &lt;dbl&gt; | regional_prob &lt;dbl&gt; | candidate_snp &lt;chr&gt; | posterior_explained_by_snp &lt;dbl&gt; | dropped_trait &lt;chr&gt; |\n",
       "|---|---|---|---|---|---|---|---|\n",
       "| 1 | 1 | 5, 9 | 0.8754 | 1 | snp917  | 0.3881 | NA |\n",
       "| 2 | 2 | 3, 8 | 0.7915 | 1 | snp2568 | 0.0418 | NA |\n",
       "| 3 | 3 | 2, 6 | 0.9426 | 1 | snp2904 | 0.1141 | NA |\n",
       "\n"
      ],
      "text/plain": [
       "  iteration traits posterior_prob regional_prob candidate_snp\n",
       "1 1         5, 9   0.8754         1             snp917       \n",
       "2 2         3, 8   0.7915         1             snp2568      \n",
       "3 3         2, 6   0.9426         1             snp2904      \n",
       "  posterior_explained_by_snp dropped_trait\n",
       "1 0.3881                     NA           \n",
       "2 0.0418                     NA           \n",
       "3 0.1141                     NA           "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "result = readRDS(\"/home/hs3393/cb_Mar/simulation_result/hyprcoloc/hyp_real_simulation_10trait/result/sample_9_real_simulation_3_ncausal_10_trait_ntr_10_hypercoloc.rds\")\n",
    "result$results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "123a9b50-968e-48aa-8d05-da26b940cae5",
   "metadata": {
    "kernel": "R",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"=====Prediction=====\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<ol>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>917</li><li>911</li><li>931</li><li>934</li><li>916</li><li>913</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>2568</li><li>2479</li><li>2605</li><li>2627</li><li>2535</li><li>2483</li><li>2485</li><li>2487</li><li>2495</li><li>2496</li><li>2497</li><li>2498</li><li>2504</li><li>2507</li><li>2510</li><li>2512</li><li>2513</li><li>2514</li><li>2515</li><li>2518</li><li>2526</li><li>2528</li><li>2529</li><li>2530</li><li>2537</li><li>2544</li><li>2545</li><li>2547</li><li>2548</li><li>2550</li><li>2552</li><li>2554</li><li>2555</li><li>2558</li><li>2560</li><li>2562</li><li>2563</li><li>2582</li><li>2590</li><li>2591</li><li>2592</li><li>2600</li><li>2603</li><li>2604</li><li>2610</li><li>2617</li><li>2618</li><li>2620</li><li>2621</li><li>2625</li><li>2630</li><li>2631</li><li>2632</li><li>2634</li><li>2635</li><li>2636</li><li>2639</li><li>2642</li><li>2644</li><li>2646</li><li>2647</li><li>2648</li><li>2650</li><li>2652</li><li>2655</li><li>2656</li><li>2657</li><li>2658</li><li>2659</li><li>2665</li><li>2667</li><li>2668</li><li>2673</li><li>2594</li><li>2539</li><li>2614</li><li>2564</li><li>2615</li><li>2649</li><li>2531</li><li>2533</li><li>2493</li><li>2576</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>2904</li><li>2892</li><li>2900</li><li>2908</li><li>2910</li><li>2911</li><li>2920</li><li>2874</li><li>2887</li><li>2859</li><li>2888</li><li>2894</li></ol>\n",
       "</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate}\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 917\n",
       "\\item 911\n",
       "\\item 931\n",
       "\\item 934\n",
       "\\item 916\n",
       "\\item 913\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 2568\n",
       "\\item 2479\n",
       "\\item 2605\n",
       "\\item 2627\n",
       "\\item 2535\n",
       "\\item 2483\n",
       "\\item 2485\n",
       "\\item 2487\n",
       "\\item 2495\n",
       "\\item 2496\n",
       "\\item 2497\n",
       "\\item 2498\n",
       "\\item 2504\n",
       "\\item 2507\n",
       "\\item 2510\n",
       "\\item 2512\n",
       "\\item 2513\n",
       "\\item 2514\n",
       "\\item 2515\n",
       "\\item 2518\n",
       "\\item 2526\n",
       "\\item 2528\n",
       "\\item 2529\n",
       "\\item 2530\n",
       "\\item 2537\n",
       "\\item 2544\n",
       "\\item 2545\n",
       "\\item 2547\n",
       "\\item 2548\n",
       "\\item 2550\n",
       "\\item 2552\n",
       "\\item 2554\n",
       "\\item 2555\n",
       "\\item 2558\n",
       "\\item 2560\n",
       "\\item 2562\n",
       "\\item 2563\n",
       "\\item 2582\n",
       "\\item 2590\n",
       "\\item 2591\n",
       "\\item 2592\n",
       "\\item 2600\n",
       "\\item 2603\n",
       "\\item 2604\n",
       "\\item 2610\n",
       "\\item 2617\n",
       "\\item 2618\n",
       "\\item 2620\n",
       "\\item 2621\n",
       "\\item 2625\n",
       "\\item 2630\n",
       "\\item 2631\n",
       "\\item 2632\n",
       "\\item 2634\n",
       "\\item 2635\n",
       "\\item 2636\n",
       "\\item 2639\n",
       "\\item 2642\n",
       "\\item 2644\n",
       "\\item 2646\n",
       "\\item 2647\n",
       "\\item 2648\n",
       "\\item 2650\n",
       "\\item 2652\n",
       "\\item 2655\n",
       "\\item 2656\n",
       "\\item 2657\n",
       "\\item 2658\n",
       "\\item 2659\n",
       "\\item 2665\n",
       "\\item 2667\n",
       "\\item 2668\n",
       "\\item 2673\n",
       "\\item 2594\n",
       "\\item 2539\n",
       "\\item 2614\n",
       "\\item 2564\n",
       "\\item 2615\n",
       "\\item 2649\n",
       "\\item 2531\n",
       "\\item 2533\n",
       "\\item 2493\n",
       "\\item 2576\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 2904\n",
       "\\item 2892\n",
       "\\item 2900\n",
       "\\item 2908\n",
       "\\item 2910\n",
       "\\item 2911\n",
       "\\item 2920\n",
       "\\item 2874\n",
       "\\item 2887\n",
       "\\item 2859\n",
       "\\item 2888\n",
       "\\item 2894\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\end{enumerate}\n"
      ],
      "text/markdown": [
       "1. 1. 917\n",
       "2. 911\n",
       "3. 931\n",
       "4. 934\n",
       "5. 916\n",
       "6. 913\n",
       "\n",
       "\n",
       "\n",
       "2. 1. 2568\n",
       "2. 2479\n",
       "3. 2605\n",
       "4. 2627\n",
       "5. 2535\n",
       "6. 2483\n",
       "7. 2485\n",
       "8. 2487\n",
       "9. 2495\n",
       "10. 2496\n",
       "11. 2497\n",
       "12. 2498\n",
       "13. 2504\n",
       "14. 2507\n",
       "15. 2510\n",
       "16. 2512\n",
       "17. 2513\n",
       "18. 2514\n",
       "19. 2515\n",
       "20. 2518\n",
       "21. 2526\n",
       "22. 2528\n",
       "23. 2529\n",
       "24. 2530\n",
       "25. 2537\n",
       "26. 2544\n",
       "27. 2545\n",
       "28. 2547\n",
       "29. 2548\n",
       "30. 2550\n",
       "31. 2552\n",
       "32. 2554\n",
       "33. 2555\n",
       "34. 2558\n",
       "35. 2560\n",
       "36. 2562\n",
       "37. 2563\n",
       "38. 2582\n",
       "39. 2590\n",
       "40. 2591\n",
       "41. 2592\n",
       "42. 2600\n",
       "43. 2603\n",
       "44. 2604\n",
       "45. 2610\n",
       "46. 2617\n",
       "47. 2618\n",
       "48. 2620\n",
       "49. 2621\n",
       "50. 2625\n",
       "51. 2630\n",
       "52. 2631\n",
       "53. 2632\n",
       "54. 2634\n",
       "55. 2635\n",
       "56. 2636\n",
       "57. 2639\n",
       "58. 2642\n",
       "59. 2644\n",
       "60. 2646\n",
       "61. 2647\n",
       "62. 2648\n",
       "63. 2650\n",
       "64. 2652\n",
       "65. 2655\n",
       "66. 2656\n",
       "67. 2657\n",
       "68. 2658\n",
       "69. 2659\n",
       "70. 2665\n",
       "71. 2667\n",
       "72. 2668\n",
       "73. 2673\n",
       "74. 2594\n",
       "75. 2539\n",
       "76. 2614\n",
       "77. 2564\n",
       "78. 2615\n",
       "79. 2649\n",
       "80. 2531\n",
       "81. 2533\n",
       "82. 2493\n",
       "83. 2576\n",
       "\n",
       "\n",
       "\n",
       "3. 1. 2904\n",
       "2. 2892\n",
       "3. 2900\n",
       "4. 2908\n",
       "5. 2910\n",
       "6. 2911\n",
       "7. 2920\n",
       "8. 2874\n",
       "9. 2887\n",
       "10. 2859\n",
       "11. 2888\n",
       "12. 2894\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[[1]]\n",
       "[1] 917 911 931 934 916 913\n",
       "\n",
       "[[2]]\n",
       " [1] 2568 2479 2605 2627 2535 2483 2485 2487 2495 2496 2497 2498 2504 2507 2510\n",
       "[16] 2512 2513 2514 2515 2518 2526 2528 2529 2530 2537 2544 2545 2547 2548 2550\n",
       "[31] 2552 2554 2555 2558 2560 2562 2563 2582 2590 2591 2592 2600 2603 2604 2610\n",
       "[46] 2617 2618 2620 2621 2625 2630 2631 2632 2634 2635 2636 2639 2642 2644 2646\n",
       "[61] 2647 2648 2650 2652 2655 2656 2657 2658 2659 2665 2667 2668 2673 2594 2539\n",
       "[76] 2614 2564 2615 2649 2531 2533 2493 2576\n",
       "\n",
       "[[3]]\n",
       " [1] 2904 2892 2900 2908 2910 2911 2920 2874 2887 2859 2888 2894\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>5</li><li>9</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>3</li><li>8</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>2</li><li>6</li></ol>\n",
       "</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate}\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 5\n",
       "\\item 9\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 3\n",
       "\\item 8\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 2\n",
       "\\item 6\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\end{enumerate}\n"
      ],
      "text/markdown": [
       "1. 1. 5\n",
       "2. 9\n",
       "\n",
       "\n",
       "\n",
       "2. 1. 3\n",
       "2. 8\n",
       "\n",
       "\n",
       "\n",
       "3. 1. 2\n",
       "2. 6\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[[1]]\n",
       "[1] 5 9\n",
       "\n",
       "[[2]]\n",
       "[1] 3 8\n",
       "\n",
       "[[3]]\n",
       "[1] 2 6\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"=====Truth=====\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>934</li><li>2485</li><li>2908</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 934\n",
       "\\item 2485\n",
       "\\item 2908\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 934\n",
       "2. 2485\n",
       "3. 2908\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1]  934 2485 2908"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<ol>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>5</li><li>6</li><li>9</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>3</li><li>8</li></ol>\n",
       "</li>\n",
       "\t<li><style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>2</li><li>6</li></ol>\n",
       "</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate}\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 5\n",
       "\\item 6\n",
       "\\item 9\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 3\n",
       "\\item 8\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\item \\begin{enumerate*}\n",
       "\\item 2\n",
       "\\item 6\n",
       "\\end{enumerate*}\n",
       "\n",
       "\\end{enumerate}\n"
      ],
      "text/markdown": [
       "1. 1. 5\n",
       "2. 6\n",
       "3. 9\n",
       "\n",
       "\n",
       "\n",
       "2. 1. 3\n",
       "2. 8\n",
       "\n",
       "\n",
       "\n",
       "3. 1. 2\n",
       "2. 6\n",
       "\n",
       "\n",
       "\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[[1]]\n",
       "[1] 5 6 9\n",
       "\n",
       "[[2]]\n",
       "[1] 3 8\n",
       "\n",
       "[[3]]\n",
       "[1] 2 6\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"=====Prediction=====\")\n",
    "# colocalization result - CoS\n",
    "result$coloc_set\n",
    "\n",
    "# colocalization result - colocalizing trait\n",
    "result$coloc_trait\n",
    "\n",
    "print(\"=====Truth=====\")\n",
    "#  true colocalizing CoS\n",
    "result$true_variant\n",
    "\n",
    "#  true colocalizing trait\n",
    "result$true_trait"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a16133c8-1d0f-4086-942b-f56ff5cad752",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Run HyprColoc without LD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dc77686-275f-4602-be54-cb927e21ca59",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[hyprcoloc_set]\n",
    "parameter: simufile = paths\n",
    "parameter: cwd = path(\"output\")\n",
    "parameter: job_size = 30\n",
    "parameter: walltime = \"50h\"\n",
    "parameter: mem = \"30G\"\n",
    "parameter: numThreads = 1\n",
    "parameter: trait = 10\n",
    "parameter: container = \"\"\n",
    "parameter: setting=\"normal\"\n",
    "input: simufile, group_by = 1\n",
    "output: f'{cwd:a}/{_input[0]:bn}_ntr_{trait}_hyprcoloc.rds'\n",
    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
    "R:  expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container \n",
    "    library(hyprcoloc)\n",
    "    source(\"/home/hs3393/cb_simulation/simulation_code/hypercoloc_set.r\")\n",
    "    file = ${_input:ar,}\n",
    "    hypercoloc_result = run_hypercoloc(file, ${trait}, setting=\"${setting}\")\n",
    "    saveRDS(hypercoloc_result, ${_output:r})\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c42da199-b853-4087-8d6b-821369689ce1",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Run Hyprcoloc with LD info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93b4d59f-c571-4f30-bec1-0653cb37a3b5",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[hyprcoloc_LD_set]\n",
    "parameter: simufile = paths\n",
    "parameter: cwd = path(\"output\")\n",
    "parameter: job_size = 30\n",
    "parameter: walltime = \"50h\"\n",
    "parameter: mem = \"30G\"\n",
    "parameter: numThreads = 1\n",
    "parameter: trait = 10\n",
    "parameter: container = \"\"\n",
    "parameter: setting=\"normal\"\n",
    "input: simufile, group_by = 1\n",
    "output: f'{cwd:a}/{_input[0]:bn}_ntr_{trait}_hypercoloc_LD.rds'\n",
    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
    "R:  expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container \n",
    "    library(hyprcoloc)\n",
    "    source(\"/home/hs3393/cb_simulation/simulation_code/hypercoloc_set_LD.r\")\n",
    "    file = ${_input:ar,}\n",
    "    hypercoloc_result = run_hypercoloc(file, ${trait}, setting=\"${setting}\")\n",
    "    saveRDS(hypercoloc_result, ${_output:r})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f752c2d9-98a8-4088-a893-617255c35c2d",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Run MOLOC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "626226ea-4421-42a0-a9c5-cf1e905597f2",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[moloc_set]\n",
    "parameter: simufile = paths\n",
    "parameter: cwd = path(\"output\")\n",
    "parameter: job_size = 30\n",
    "parameter: walltime = \"50h\"\n",
    "parameter: mem = \"40G\"\n",
    "parameter: numThreads = 1\n",
    "parameter: trait = 50\n",
    "parameter: container = \"\"\n",
    "input: simufile, group_by = 1\n",
    "output: f'{cwd:a}/{_input[0]:bn}_ntr_{trait}_moloc_set.rds'\n",
    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
    "R:  expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container \n",
    "    library(moloc)\n",
    "    source(\"/home/hs3393/cb_simulation/simulation_code/moloc_set.r\")\n",
    "    files = c(${_input:ar,})\n",
    "    moloc_result = run_moloc(files, ${trait})\n",
    "    saveRDS(moloc_result, ${_output:r})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "333a696d-ed5b-45fe-869e-24195d77d98c",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Run COLOC V5 (SuSiE coloc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "830a4261-c51d-4e15-9641-520ec45532a6",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "[susie_coloc_set]\n",
    "parameter: simufile = paths\n",
    "parameter: cwd = path(\"output\")\n",
    "parameter: job_size = 300\n",
    "parameter: walltime = \"50h\"\n",
    "parameter: mem = \"40G\"\n",
    "parameter: numThreads = 1\n",
    "parameter: trait = 2\n",
    "parameter: container = \"\"\n",
    "input: simufile, group_by = 1\n",
    "output: f'{cwd:a}/{_input[0]:bn}_ntr_{trait}_susie_coloc_result.rds'\n",
    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'\n",
    "R:  expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container \n",
    "    source(\"/home/hs3393/cb_simulation/simulation_code/susie_coloc.r\")\n",
    "    file = ${_input:ar,}\n",
    "    susie_coloc_result = run_susie_coloc(file)\n",
    "    saveRDS(susie_coloc_result, ${_output:r})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6c9527f-3b79-48a9-b77f-89a42e6233e1",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "## Bash files\n",
    "\n",
    "### Hyprcoloc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b7d8745-09ca-4886-9302-68dd0589c838",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "## 2 trait\n",
    "\n",
    "data_dir=\"/home/hs3393/cb_Mar/simulation_data/real_simulation_2trait/\"\n",
    "job=\"hyp_real_simulation_2trait\"\n",
    "work_dir=\"/home/hs3393/cb_Mar/simulation_result/hyprcoloc/\"\n",
    "#!/bin/bash\n",
    "\n",
    "mkdir -p ${work_dir}/${job}/code\n",
    "mkdir -p ${work_dir}/${job}/log\n",
    "mkdir -p ${work_dir}/${job}/result\n",
    "\n",
    "cd ${work_dir}/${job}/code\n",
    "\n",
    "cat << 'EOF' > base_script\n",
    "#!/bin/bash -l\n",
    "# NOTE the -l flag!\n",
    "#\n",
    "#SBATCH -t 80:00:00\n",
    "#SBATCH --mem=20000\n",
    "#SBATCH -J JOB\n",
    "#SBATCH -o WORK_DIR/JOB/log/JOB.%j.out\n",
    "#SBATCH -e WORK_DIR/JOB/log/JOB.%j.err\n",
    "\n",
    "source /home/hs3393/mamba_activate.sh\n",
    "module load Singularity\n",
    "\n",
    "cd DATA_DIR\n",
    "sos run /home/hs3393/cb_Mar/simulation_code/3.Other_Methods.ipynb hyprcoloc_set \\\n",
    "    --simufile $(find -type f -name '*_NCAUSAL_ncausal_*.rds') \\\n",
    "    --mem 20G --trait 2 \\\n",
    "    --cwd WORK_DIR/JOB/result\n",
    "EOF\n",
    "\n",
    "\n",
    "base_script=\"base_script\"\n",
    "for ncausal in 1 2 3; do\n",
    "    output_script=\"ncausal_${ncausal}.sh\"\n",
    "    cat ${base_script}| sed \"s|WORK_DIR|${work_dir}|g\" |sed \"s|NCAUSAL|${ncausal}|g\" | sed \"s|JOB|${job}|g\" | sed \"s|DATA_DIR|${data_dir}|g\"  > ${output_script}\n",
    "    sbatch ${output_script}\n",
    "done\n",
    "\n",
    "## 5 trait\n",
    "\n",
    "data_dir=\"/home/hs3393/cb_Mar/simulation_data/real_simulation_5trait/\"\n",
    "job=\"hyp_real_simulation_5trait\"\n",
    "work_dir=\"/home/hs3393/cb_Mar/simulation_result/hyprcoloc/\"\n",
    "#!/bin/bash\n",
    "\n",
    "mkdir -p ${work_dir}/${job}/code\n",
    "mkdir -p ${work_dir}/${job}/log\n",
    "mkdir -p ${work_dir}/${job}/result\n",
    "\n",
    "cd ${work_dir}/${job}/code\n",
    "\n",
    "cat << 'EOF' > base_script\n",
    "#!/bin/bash -l\n",
    "# NOTE the -l flag!\n",
    "#\n",
    "#SBATCH -t 80:00:00\n",
    "#SBATCH --mem=20000\n",
    "#SBATCH -J JOB\n",
    "#SBATCH -o WORK_DIR/JOB/log/JOB.%j.out\n",
    "#SBATCH -e WORK_DIR/JOB/log/JOB.%j.err\n",
    "\n",
    "source /home/hs3393/mamba_activate.sh\n",
    "module load Singularity\n",
    "\n",
    "cd DATA_DIR\n",
    "sos run /home/hs3393/cb_Mar/simulation_code/3.Other_Methods.ipynb hyprcoloc_set \\\n",
    "    --simufile $(find -type f -name '*_NCAUSAL_ncausal_*.rds') \\\n",
    "    --mem 20G --trait 5 \\\n",
    "    --cwd WORK_DIR/JOB/result\n",
    "EOF\n",
    "\n",
    "\n",
    "base_script=\"base_script\"\n",
    "for ncausal in 1 2 3 4; do\n",
    "    output_script=\"ncausal_${ncausal}.sh\"\n",
    "    cat ${base_script}| sed \"s|WORK_DIR|${work_dir}|g\" |sed \"s|NCAUSAL|${ncausal}|g\" | sed \"s|JOB|${job}|g\" | sed \"s|DATA_DIR|${data_dir}|g\"  > ${output_script}\n",
    "    sbatch ${output_script}\n",
    "done\n",
    "\n",
    "## 10 trait\n",
    "\n",
    "data_dir=\"/home/hs3393/cb_Mar/simulation_data/real_simulation_10trait/\"\n",
    "job=\"hyp_real_simulation_10trait\"\n",
    "work_dir=\"/home/hs3393/cb_Mar/simulation_result/hyprcoloc/\"\n",
    "#!/bin/bash\n",
    "\n",
    "mkdir -p ${work_dir}/${job}/code\n",
    "mkdir -p ${work_dir}/${job}/log\n",
    "mkdir -p ${work_dir}/${job}/result\n",
    "\n",
    "cd ${work_dir}/${job}/code\n",
    "\n",
    "cat << 'EOF' > base_script\n",
    "#!/bin/bash -l\n",
    "# NOTE the -l flag!\n",
    "#\n",
    "#SBATCH -t 80:00:00\n",
    "#SBATCH --mem=30000\n",
    "#SBATCH -J JOB\n",
    "#SBATCH -o WORK_DIR/JOB/log/JOB.%j.out\n",
    "#SBATCH -e WORK_DIR/JOB/log/JOB.%j.err\n",
    "\n",
    "source /home/hs3393/mamba_activate.sh\n",
    "module load Singularity\n",
    "\n",
    "cd DATA_DIR\n",
    "sos run /home/hs3393/cb_Mar/simulation_code/3.Other_Methods.ipynb hyprcoloc_set \\\n",
    "    --simufile $(find -type f -name '*_NCAUSAL_ncausal_*.rds') \\\n",
    "    --mem 30G --trait 10 \\\n",
    "    --cwd WORK_DIR/JOB/result\n",
    "EOF\n",
    "\n",
    "\n",
    "base_script=\"base_script\"\n",
    "for ncausal in 1 2 3 4; do\n",
    "    output_script=\"ncausal_${ncausal}.sh\"\n",
    "    cat ${base_script}| sed \"s|WORK_DIR|${work_dir}|g\" |sed \"s|NCAUSAL|${ncausal}|g\" | sed \"s|JOB|${job}|g\" | sed \"s|DATA_DIR|${data_dir}|g\"  > ${output_script}\n",
    "    sbatch ${output_script}\n",
    "done\n",
    "\n",
    "## 20 trait\n",
    "\n",
    "data_dir=\"/home/hs3393/cb_Mar/simulation_data/real_simulation_20trait/\"\n",
    "job=\"hyp_real_simulation_20trait\"\n",
    "work_dir=\"/home/hs3393/cb_Mar/simulation_result/hyprcoloc/\"\n",
    "#!/bin/bash\n",
    "\n",
    "mkdir -p ${work_dir}/${job}/code\n",
    "mkdir -p ${work_dir}/${job}/log\n",
    "mkdir -p ${work_dir}/${job}/result\n",
    "\n",
    "cd ${work_dir}/${job}/code\n",
    "\n",
    "cat << 'EOF' > base_script\n",
    "#!/bin/bash -l\n",
    "# NOTE the -l flag!\n",
    "#\n",
    "#SBATCH -t 80:00:00\n",
    "#SBATCH --mem=40000\n",
    "#SBATCH -J JOB\n",
    "#SBATCH -o WORK_DIR/JOB/log/JOB.%j.out\n",
    "#SBATCH -e WORK_DIR/JOB/log/JOB.%j.err\n",
    "\n",
    "source /home/hs3393/mamba_activate.sh\n",
    "module load Singularity\n",
    "\n",
    "cd DATA_DIR\n",
    "sos run /home/hs3393/cb_Mar/simulation_code/3.Other_Methods.ipynb hyprcoloc_set \\\n",
    "    --simufile $(find -type f -name '*_NCAUSAL_ncausal_*.rds') \\\n",
    "    --mem 40G --trait 20 \\\n",
    "    --cwd WORK_DIR/JOB/result\n",
    "EOF\n",
    "\n",
    "base_script=\"base_script\"\n",
    "for ncausal in 1 2 3 4; do\n",
    "    output_script=\"ncausal_${ncausal}.sh\"\n",
    "    cat ${base_script}| sed \"s|WORK_DIR|${work_dir}|g\" |sed \"s|NCAUSAL|${ncausal}|g\" | sed \"s|JOB|${job}|g\" | sed \"s|DATA_DIR|${data_dir}|g\"  > ${output_script}\n",
    "    sbatch ${output_script}\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82e188c0-027f-409b-951b-066c4455af54",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "### MOLOC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3abad90a-f29c-4a08-ab59-e4768564d43c",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "data_dir=\"/home/hs3393/cb_Mar/simulation_data/real_simulation_2trait/\"\n",
    "job=\"moloc_real_simulation_2trait\"\n",
    "work_dir=\"/home/hs3393/cb_Mar/simulation_result/moloc\"\n",
    "#!/bin/bash\n",
    "\n",
    "mkdir -p ${work_dir}/${job}/code\n",
    "mkdir -p ${work_dir}/${job}/log\n",
    "mkdir -p ${work_dir}/${job}/result\n",
    "\n",
    "cd ${work_dir}/${job}/code\n",
    "\n",
    "cat << 'EOF' > base_script\n",
    "#!/bin/bash -l\n",
    "# NOTE the -l flag!\n",
    "#\n",
    "#SBATCH -t 80:00:00\n",
    "#SBATCH --mem=15000\n",
    "#SBATCH -J JOB\n",
    "#SBATCH -o WORK_DIR/JOB/log/JOB.%j.out\n",
    "#SBATCH -e WORK_DIR/JOB/log/JOB.%j.err\n",
    "\n",
    "source /home/hs3393/mamba_activate.sh\n",
    "module load Singularity\n",
    "\n",
    "cd DATA_DIR\n",
    "sos run /home/hs3393/cb_Mar/simulation_code/3.Other_Methods.ipynb moloc_set \\\n",
    "    --simufile $(find -type f -name '*_NCAUSAL_ncausal_*.rds') \\\n",
    "    --mem 15G --trait 2 \\\n",
    "    --cwd WORK_DIR/JOB/result\n",
    "EOF\n",
    "\n",
    "\n",
    "base_script=\"base_script\"\n",
    "for ncausal in 1 2 3; do\n",
    "    output_script=\"ncausal_${ncausal}.sh\"\n",
    "    cat ${base_script}| sed \"s|WORK_DIR|${work_dir}|g\" |sed \"s|NCAUSAL|${ncausal}|g\" | sed \"s|JOB|${job}|g\" | sed \"s|DATA_DIR|${data_dir}|g\"  > ${output_script}\n",
    "    sbatch ${output_script}\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08145458-2cb0-4a07-8d26-2cfb4bedb326",
   "metadata": {
    "kernel": "SoS",
    "tags": []
   },
   "source": [
    "### COLOC V5 (SuSiE-COLOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d2524ce-b77b-409f-9cf7-7066c855d228",
   "metadata": {
    "kernel": "SoS"
   },
   "outputs": [],
   "source": [
    "data_dir=\"/home/hs3393/cb_Mar/simulation_data/real_simulation_2trait/\"\n",
    "job=\"susie_coloc_real_simulation_2trait\"\n",
    "work_dir=\"/home/hs3393/cb_Mar/simulation_result/susie_coloc\"\n",
    "#!/bin/bash\n",
    "\n",
    "mkdir -p ${work_dir}/${job}/code\n",
    "mkdir -p ${work_dir}/${job}/log\n",
    "mkdir -p ${work_dir}/${job}/result\n",
    "\n",
    "cd ${work_dir}/${job}/code\n",
    "\n",
    "cat << 'EOF' > base_script\n",
    "#!/bin/bash -l\n",
    "# NOTE the -l flag!\n",
    "#\n",
    "#SBATCH -t 80:00:00\n",
    "#SBATCH --mem=40000\n",
    "#SBATCH -J JOB\n",
    "#SBATCH -o WORK_DIR/JOB/log/JOB.%j.out\n",
    "#SBATCH -e WORK_DIR/JOB/log/JOB.%j.err\n",
    "\n",
    "source /home/hs3393/mamba_activate.sh\n",
    "module load Singularity\n",
    "\n",
    "cd DATA_DIR\n",
    "sos run /home/hs3393/cb_Mar/simulation_code/3.Other_Methods.ipynb susie_coloc_set \\\n",
    "    --simufile $(find -type f -name '*_NCAUSAL_ncausal_*.rds') \\\n",
    "    --mem 40G --trait 2 \\\n",
    "    --cwd WORK_DIR/JOB/result\n",
    "EOF\n",
    "\n",
    "\n",
    "base_script=\"base_script\"\n",
    "for ncausal in 1 2 3; do\n",
    "    output_script=\"ncausal_${ncausal}.sh\"\n",
    "    cat ${base_script}| sed \"s|WORK_DIR|${work_dir}|g\" |sed \"s|NCAUSAL|${ncausal}|g\" | sed \"s|JOB|${job}|g\" | sed \"s|DATA_DIR|${data_dir}|g\"  > ${output_script}\n",
    "    sbatch ${output_script}\n",
    "done"
   ]
  }
 ],
 "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": {
   "kernels": [
    [
     "R",
     "ir",
     "R",
     "",
     ""
    ],
    [
     "SoS",
     "sos",
     "",
     "",
     "sos"
    ]
   ],
   "version": "0.24.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}