Merge pull request #1095 from Chunmingl/main

update for parameters and corner cases

code/pecotmr_integration/twas.ipynb

@@ -174,9 +174,9 @@
     "   --gwas_meta_data /mnt/vast/hpc/csg/cl4215/mrmash/workflow/GWAS/gwas_meta.tsv \\\n",
     "   --ld_meta_data /mnt/vast/hpc/csg/data_public/20240409_ADSP_LD_matrix/ld_meta_file.tsv \\\n",
     "   --regions /mnt/vast/hpc/csg/cl4215/mrmash/workflow/pipeline_data/picalm_EUR_LD_blocks.bed \\\n",
-    "   --xqtl_meta_data /mnt/vast/hpc/csg/cl4215/mrmash/workflow/pipeline_data/picalm_region_gene_meta_data_test.tsv \\\n",
+    "   --xqtl_meta_data /mnt/vast/hpc/csg/cl4215/mrmash/workflow/pipeline_data/twas_157_genes_table_2024.tsv \\\n",
     "   --xqtl_type_table /mnt/vast/hpc/csg/cl4215/mrmash/workflow/pipeline_data/data_type_table.txt \\\n",
-    "   --p_value_cutoff 0.05 --rsq_threshold 0.01\n",
+    "   --rsq_pval_cutoff 0.05 --rsq_cutoff 0.01\n",
     "```"
    ]
   },
@@ -433,11 +433,13 @@
     "depends: sos_variable(\"filtered_regional_xqtl_files\")\n",
     "parameter: coverage = \"cs_coverage_0.95\"\n",
     "# Threashold for rsq and pvalue for imputability determination for a model \n",
-    "parameter: rsq_threshold = 0.01\n",
-    "parameter: rsq_pval_threshold = 0.05\n",
+    "parameter: rsq_cutoff = 0.01\n",
+    "parameter: rsq_pval_cutoff = 0.05\n",
     "parameter: mr_pval_cutoff = 0.05\n",
     "parameter: save_ctwas_data = True\n",
     "parameter: save_mr_result = True\n",
+    "parameter: rsq_option=\"adj_rsq\"\n",
+    "parameter: rsq_pval_option=[\"adj_rsq_pval\", \"pval\"]\n",
     "\n",
     "input: filtered_regional_xqtl_files, group_by = lambda x: group_by_region(x, filtered_regional_xqtl_files), group_with = \"filtered_region_info\"\n",
     "output_files = [f'{cwd:a}/{step_name}/{name}.{_filtered_region_info[3]}.twas.tsv.gz']\n",
@@ -485,14 +487,19 @@
     "                                     \"${ld_meta_data}\", \n",
     "                                     \"${gwas_meta_data}\",  \n",
     "                                     region_block=\"${_filtered_region_info[3]}\",\n",
-    "                                     rsq_threshold = ${rsq_threshold}, \n",
-    "                                     rsq_pval_threshold = ${p_value_cutoff}, \n",
+    "                                     rsq_cutoff = ${rsq_cutoff}, \n",
+    "                                     rsq_option = \"${rsq_option}\", \n",
+    "                                     rsq_pval_cutoff = ${rsq_pval_cutoff}, \n",
+    "                                     rsq_pval_option=c(${\", \".join([f'\"{x}\"' for x in rsq_pval_option])}), \n",
     "                                     mr_pval_cutoff = ${mr_pval_cutoff},\n",
-    "                                     mr_coverage_column = ${coverage}, \n",
+    "                                     mr_coverage_column = \"${coverage}\", \n",
     "                                     output_twas_data = ${\"TRUE\" if save_ctwas_data else \"FALSE\"})\n",
-    "    twas_results_db$twas_result <- merge( twas_results_db$twas_result, xqtl_meta_df[, c(\"region_id\",\"TSS\",\"start\",\"end\")], by.x=\"molecular_id\", by.y=\"region_id\")\n",
-    "    fwrite(twas_results_db$twas_resul[, c(2,1,14:16,3:13)], file = ${_output[0]:r}, sep = \"\\t\", compress = \"gzip\")\n",
-    "\n",
+    "    if(!is.null(twas_results_db)){\n",
+    "        twas_results_db$twas_result <- merge( twas_results_db$twas_result, xqtl_meta_df[, c(\"region_id\",\"TSS\",\"start\",\"end\")], by.x=\"molecular_id\", by.y=\"region_id\")\n",
+    "        fwrite(twas_results_db$twas_result[, c(2,1,14:16,3:13)], file = ${_output[0]:r}, sep = \"\\t\", compress = \"gzip\")\n",
+    "    } else {\n",
+    "        fwrite(data.frame(), file = ${_output[0]:r}, sep = \"\\t\", compress = \"gzip\")\n",
+    "    }\n",
     "    # Step 3: reformat for follow up cTWAS analysis\n",
     "    if (${\"TRUE\" if save_ctwas_data else \"FALSE\"}) {\n",
     "        saveRDS(twas_results_db$twas_data, \"${_output[0]:nnn}.twas_data.rds\", compress='xz')\n",