Add a dynamic way to increase numer of xQTL to map in fine-mapping

cumc · Oct 24, 2023 · 799f98a · 799f98a
1 parent 3514da4
commit 799f98a
Showing 1 changed file with 21 additions and 8 deletions.
diff --git a/code/fine_mapping/SuSiE/SuSiE.ipynb b/code/fine_mapping/SuSiE/SuSiE.ipynb
@@ -330,7 +330,8 @@
    "outputs": [],
    "source": [
     "[susie_1]\n",
-    "parameter: max_L = 20\n",
+    "parameter: init_L = 12\n",
+    "parameter: max_L = 30\n",
     "# remove a variant if it has more than imiss missing individual level data\n",
     "parameter: imiss = 1\n",
     "# MAF cutoff\n",
@@ -366,21 +367,33 @@
     "                                          keep_indel = ${\"TRUE\" if indel else \"FALSE\"})\n",
     "   \n",
     "    # Fine-mapping with SuSiE\n",
-    "    library(susieR)  \n",
+    "    library(susieR)\n",
+    "    l_step = 5\n",
     "    fitted = list()\n",
     "    for (r in 1:length(fdat$residual_Y_scaled)) { ## Cant have a universal way to specify names due to the accomodation of missingness, use index instead\n",
-    "        st = proc.time()\n",
-    "        fitted[[r]] <- susie(fdat$residual_X_scaled[[r]],\n",
+    "        L = ${init_L}\n",
+    "        # Dynamically increase L\n",
+    "        while (TRUE) {\n",
+    "            st = proc.time()\n",
+    "            fitted[[r]] <- susie(fdat$residual_X_scaled[[r]],\n",
     "                             fdat$residual_Y_scaled[[r]],\n",
-    "                             L=${max_L},\n",
+    "                             L=L,\n",
     "                             max_iter=500,\n",
     "                             estimate_residual_variance=TRUE,\n",
     "                             estimate_prior_variance=TRUE,\n",
     "                             refine=TRUE,\n",
     "                             compute_univariate_zscore=FALSE,\n",
     "                             coverage=${coverage})\n",
-    "        fitted[[r]]$analysis_time <- proc.time() - st\n",
-    "        fitted[[r]] <- post_process_susie(fitted[[r]], fdat, r, signal_cutoff = ${pip_cutoff})\n",
+    "            fitted[[r]]$analysis_time <- proc.time() - st\n",
+    "            fitted[[r]] <- post_process_susie(fitted[[r]], fdat, r, signal_cutoff = ${pip_cutoff})\n",
+    "            if (!is.null(fitted[[r]]$sets)) {\n",
+    "                if (length(fitted[[r]]$sets)>=L && L<=${max_L}) {\n",
+    "                  L = L + l_step\n",
+    "                } else {\n",
+    "                  break\n",
+    "                }\n",
+    "            }\n",
+    "        }\n",
     "    }\n",
     "    names(fitted) <- names(fdat$residual_Y_scaled)\n",
     "    saveRDS(fitted, ${_output:ar}, compress='xz')"
@@ -674,7 +687,7 @@
      ""
     ]
    ],
-   "version": "0.22.4"
+   "version": "0.20.1"
   }
  },
  "nbformat": 4,