diff --git a/code/mnm_analysis/mnm_methods/colocboost.ipynb b/code/mnm_analysis/mnm_methods/colocboost.ipynb index 4e769c005..c8e61d129 100644 --- a/code/mnm_analysis/mnm_methods/colocboost.ipynb +++ b/code/mnm_analysis/mnm_methods/colocboost.ipynb @@ -291,7 +291,9 @@ }, { "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ "# Multi-ancestry: per-study LD reference specified in gwas_meta_multi.tsv\n", "sos run pipeline/colocboost.ipynb colocboost \\\n", @@ -306,9 +308,7 @@ " --separate-gwas --xqtl-coloc \\\n", " --region-name ENSG00000239945 \\\n", " --phenotype-names trait_A\n" - ], - "outputs": [], - "execution_count": null + ] }, { "cell_type": "markdown", @@ -1182,7 +1182,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": { "kernel": "SoS", "tags": [] @@ -1350,21 +1350,21 @@ " \n", " # - analysis parameters\n", " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", - " res <- colocboost_analysis_pipeline(region_data = fdat, \n", - " event_filters = event_filters,\n", - " # - analysis\n", - " xqtl_coloc = xqtl_coloc,\n", - " joint_gwas = joint_gwas,\n", - " separate_gwas = separate_gwas,\n", - " # - individual QC\n", - " maf_cutoff = maf_cutoff, \n", - " pip_cutoff_to_skip_ind = pip_cutoff_to_skip_ind,\n", - " # - sumstat QC\n", - " pip_cutoff_to_skip_sumstat = pip_cutoff_to_skip_sumstat,\n", - " keep_indel = keep_indel,\n", - " qc_method = qc_method,\n", - " impute = impute, \n", - " impute_opts = impute_opts)\n", + " res <- colocboost_pipeline(region_data = fdat, \n", + " event_filters = event_filters,\n", + " # - analysis\n", + " xqtl_coloc = xqtl_coloc,\n", + " joint_gwas = joint_gwas,\n", + " separate_gwas = separate_gwas,\n", + " # - individual QC\n", + " maf_cutoff = maf_cutoff, \n", + " pip_cutoff_to_skip_ind = pip_cutoff_to_skip_ind,\n", + " # - sumstat QC\n", + " pip_cutoff_to_skip_sumstat = pip_cutoff_to_skip_sumstat,\n", + " keep_indel = keep_indel,\n", + " qc_method = qc_method,\n", + " impute = impute, \n", + " impute_opts = impute_opts)\n", "\n", " # Reorganize outputs\n", " computing_time <- c(computing_time, res$computing_time)\n", diff --git a/code/mnm_analysis/mnm_methods/mnm_regression.ipynb b/code/mnm_analysis/mnm_methods/mnm_regression.ipynb index c6a5b1d73..d6c6e1197 100644 --- a/code/mnm_analysis/mnm_methods/mnm_regression.ipynb +++ b/code/mnm_analysis/mnm_methods/mnm_regression.ipynb @@ -944,7 +944,7 @@ " region_name = \"${_meta_info[2]}\"\n", " }\n", "\n", - " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", + " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), analysis_region = \"${_meta_info[1]}\", grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", "\n", " finemapping_result = list()\n", " preset_variants_result = list()\n", @@ -979,7 +979,8 @@ " X_scalar=fdat$residual_X_scalar[[r]], \n", " Y_scalar=if (fdat$residual_Y_scalar[[r]] == 1) 1 else fdat$residual_Y_scalar[[r]][,i,drop=FALSE],\n", " X_variance=fdat$X_variance[[r]],\n", - " other_quantities = list(dropped_samples = dropped_samples),\n", + " region = region_info$analysis_region,\n", + " other_quantities = list(dropped_samples = dropped_samples, condition_id = names(fdat$residual_Y)[r]),\n", " # filters\n", " imiss_cutoff = ${imiss},\n", " maf_cutoff = NULL,\n", @@ -1014,7 +1015,8 @@ " X_scalar=fdat$residual_X_scalar[[r]], \n", " Y_scalar=if (fdat$residual_Y_scalar[[r]] == 1) 1 else fdat$residual_Y_scalar[[r]][,i,drop=FALSE],\n", " X_variance=fdat$X_variance[[r]],\n", - " other_quantities = list(dropped_samples = dropped_samples),\n", + " region = region_info$analysis_region,\n", + " other_quantities = list(dropped_samples = dropped_samples, condition_id = names(fdat$residual_Y)[r]),\n", " # filters\n", " imiss_cutoff = ${imiss},\n", " maf_cutoff = ${min_twas_maf},\n", @@ -1232,7 +1234,7 @@ "\n", " context_regions <- list(${\",\".join([(\"c('\"+x+\"')\") if isinstance(x, str) else (\"c\"+ str(x)) for x in _meta_info[3]])})\n", " names(context_regions) <- colnames(fdat$residual_Y)\n", - " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", + " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), analysis_region = \"${_meta_info[1]}\", grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", "\n", " dd_prior <- if (${mixture_prior:r} == '.' || ${mixture_prior:r} == '') NULL else readRDS(${mixture_prior:r})\n", " dd_prior_cv <- if (${mixture_prior_cv:r} == '.' || ${mixture_prior_cv:r} == '') NULL else readRDS(${mixture_prior_cv:r})\n", @@ -1244,6 +1246,7 @@ " Y = fdat$residual_Y,\n", " maf = fdat$maf,\n", " X_variance = fdat$X_variance,\n", + " region = region_info$analysis_region,\n", " other_quantities = list(dropped_samples = fdat$dropped_samples),\n", " # filters\n", " imiss_cutoff = ${imiss}, \n", @@ -1414,14 +1417,14 @@ " gene_id <- names(combined_list[[\"extracted_regional_window_combined_susie_result\"]][[condition_name]])[i]\n", " \n", " # Extract top_loci variants\n", - " variants1 <- unique(susie_result_1[[\"top_loci\"]][[\"variant_id\"]])\n", + " variants1 <- unique(susie_result_1[[\"top_loci\"]][[\"variant\"]])\n", " if (length(variants1) > 0) {\n", " different_in_all_genes <- TRUE\n", " \n", " for (j in 1:length(combined_list[[\"extracted_regional_window_combined_susie_result\"]][[condition_name]])) {\n", " if (i != j) {\n", " susie_result_2 <- combined_list[[\"extracted_regional_window_combined_susie_result\"]][[condition_name]][[j]]\n", - " variants2 <- unique(susie_result_2[[\"top_loci\"]][[\"variant_id\"]])\n", + " variants2 <- unique(susie_result_2[[\"top_loci\"]][[\"variant\"]])\n", " \n", " # Compare variants\n", " if (length(variants2) > 0) {\n", @@ -1515,7 +1518,7 @@ " region_name = \"${_meta_info[2]}\"\n", " }\n", "\n", - " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", + " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), analysis_region = \"${_meta_info[1]}\", grange=parse_region(\"${_meta_info[1]}\"), region_name=region_name)\n", " \n", " dd_prior <- NULL\n", " if (${\"TRUE\" if data_driven_prior else \"FALSE\"}) {\n", @@ -1532,6 +1535,7 @@ " Y=fdat$residual_Y[,filtered_event_names],\n", " maf=fdat$maf,\n", " X_variance = fdat$X_variance,\n", + " region = region_info$analysis_region,\n", " other_quantities = list(dropped_samples = fdat$dropped_samples),\n", " # filters\n", " imiss_cutoff = ${imiss}, \n", @@ -1664,6 +1668,7 @@ " }\n", " \n", " fitted = setNames(replicate(length(fdat$residual_Y), list(), simplify = FALSE), names(fdat$residual_Y))\n", + " region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), analysis_region = \"${_meta_info[1]}\", grange=parse_region(\"${_meta_info[1]}\"), region_name=\"${_meta_info[2]}\")\n", " for (r in 1:length(fitted)) {\n", " st = proc.time()\n", " fitted[[r]] = list()\n", @@ -1684,6 +1689,7 @@ " data_x = fdat$residual_X[[r]], data_y = top_pc_data[,i],\n", " X_scalar = fdat$residual_X_scalar[[r]], y_scalar = 1, maf = fdat$maf[[r]],\n", " coverage = ${coverage[0]}, secondary_coverage = c(${\",\".join([str(x) for x in coverage[1:]])}), signal_cutoff = ${pip_cutoff},\n", + " region = region_info$analysis_region,\n", " other_quantities = list(dropped_samples = list(X=fdat$dropped_sample$dropped_samples_X[[r]], \n", " y=fdat$dropped_sample$dropped_samples_Y[[r]], \n", " covar=fdat$dropped_sample$dropped_samples_covar[[r]]))),\n", @@ -1722,11 +1728,12 @@ " data_x = fdat$residual_X[[r]], data_y = NULL,\n", " X_scalar = fdat$residual_X_scalar[[r]], y_scalar = 1, maf = fdat$maf[[r]],\n", " coverage = ${coverage[0]}, secondary_coverage = c(${\",\".join([str(x) for x in coverage[1:]])}), signal_cutoff = ${pip_cutoff},\n", + " region = region_info$analysis_region,\n", " other_quantities = list(dropped_samples = list(X=fdat$dropped_sample$dropped_samples_X[[r]], y=fdat$dropped_sample$dropped_samples_Y[[r]], \n", " covar=fdat$dropped_sample$dropped_samples_covar[[r]]))),\n", " primary_method = \"fsusie\")\n", " fitted[[r]]$total_time_elapsed = proc.time() - st\n", - " fitted[[r]]$region_info = list(region_coord=parse_region(\"${_meta_info[0]}\"), grange=parse_region(\"${_meta_info[1]}\"), region_name=\"${_meta_info[2]}\")\n", + " fitted[[r]]$region_info = region_info\n", " # original data no longer relevant, set to NA to release memory\n", " fdat$residual_X[[r]] <- NA\n", " fdat$residual_Y[[r]] <- NA\n", @@ -1827,4 +1834,4 @@ }, "nbformat": 4, "nbformat_minor": 4 -} +} \ No newline at end of file