Skip to content

twas, susie_enloc and ctwas update #1381

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
al4225 merged 5 commits into from
Jun 25, 2026
Merged

twas, susie_enloc and ctwas update #1381

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0b95013
twas update
danielnachun Jun 24, 2026
a5f7cb2
twas update
danielnachun Jun 24, 2026
2e1d6fd
enloc updates
danielnachun Jun 24, 2026
c7df484
add ctwas fixes
danielnachun Jun 24, 2026
41268ae
add last ctwas bits
danielnachun Jun 24, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
139 changes: 52 additions & 87 deletions code/SoS/pecotmr_integration/SuSiE_enloc.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -561,19 +561,31 @@
"input: regional_data[\"data\"], group_by = lambda x: group_by_region(x, regional_data[\"data\"]), group_with = \"meta\"\n",
"output: f'{cwd:a}/{name}.{_meta[0]}.enrichment.rds'\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
"R: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container\n",
" library(tidyverse)\n",
" library(pecotmr)\n",
" # RDS files for GWAS data\n",
" gwas_finemapped_data = c(${paths([x for x in _meta[1:len(_meta)]]):r,})\n",
" # RDS files for xQTL data\n",
" xqtl_finemapped_data = c(${paths([x for x in _input]):r,})\n",
" enrich_result = xqtl_enrichment_wrapper(gwas_files = gwas_finemapped_data, xqtl_files = xqtl_finemapped_data, \n",
" xqtl_finemapping_obj = c(\"${_meta[0]}\",${\",\".join(['\"%s\"' % x for x in xqtl_finemapping_obj]) if len(xqtl_finemapping_obj) != 0 else \"NULL\"}), \n",
" xqtl_varname_obj = c(\"${_meta[0]}\",${\",\".join(['\"%s\"' % x for x in xqtl_varname_obj]) if len(xqtl_varname_obj) != 0 else \"NULL\"}), \n",
" gwas_finemapping_obj = c(${\",\".join(['\"%s\"' % x for x in gwas_finemapping_obj]) if len(gwas_finemapping_obj) != 0 else \"NULL\"}), \n",
" gwas_varname_obj = c(${\",\".join(['\"%s\"' % x for x in gwas_varname_obj]) if len(gwas_varname_obj) != 0 else \"NULL\"}))\n",
" saveRDS(enrich_result, ${_output:ar})"
"bash: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container\n",
" set -e\n",
" # Per-region xQTL files for this group, and the GWAS rds entries in _meta[1:].\n",
" # Convert each side of this region to a pecotmr S4 collection, then run\n",
" # qtlEnrichmentPipeline on the pair.\n",
" Rscript code/script/pecotmr_integration/legacy_enloc_finemap_convert.R \\\n",
" --mode qtl \\\n",
" --rds-files ${\",\".join([str(x) for x in _input])} \\\n",
" --finemapping-obj '${\" \".join(xqtl_finemapping_obj)}' \\\n",
" --varname-obj '${\" \".join(xqtl_varname_obj)}' \\\n",
" --study ${name} \\\n",
" --output ${_output:n}.qtl.rds\n",
"\n",
" Rscript code/script/pecotmr_integration/legacy_enloc_finemap_convert.R \\\n",
" --mode gwas \\\n",
" --rds-files ${\",\".join([str(x) for x in _meta[1:] if str(x).endswith(\".rds\")])} \\\n",
" --finemapping-obj '${\" \".join(gwas_finemapping_obj)}' \\\n",
" --varname-obj '${\" \".join(gwas_varname_obj)}' \\\n",
" --output ${_output:n}.gwas.rds\n",
"\n",
" Rscript code/script/pecotmr_integration/qtl_enrichment.R \\\n",
" --qtl-fine-mapping ${_output:n}.qtl.rds \\\n",
" --gwas-fine-mapping ${_output:n}.gwas.rds \\\n",
" --ncore ${numThreads} \\\n",
" --output ${_output}\n"
]
},
{
Expand All @@ -600,79 +612,32 @@
"# output: f'{cwd:a}/{step_name[:-2]}/{name}.{_meta[0]}.coloc.rds'\n",
"output: f'{cwd:a}/{step_name}/{name}.{_meta[0]}.coloc.rds'\n",
"task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'\n",
"R: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container\n",
" library(tidyverse)\n",
" library(pecotmr)\n",
" library(coloc) \n",
" # RDS files for xQTL data\n",
" # RDS files for xQTL data\n",
" xqtl_finemapped_datas = c(${paths([x for x in _input]):r,})\n",
" coloc_res <- list()\n",
" # are we doing something like this? that means results are by each context, and each one has \n",
" for(xqtl_finemapped_data in xqtl_finemapped_datas){\n",
" qtl_dat <- readRDS(xqtl_finemapped_data)\n",
" gene = names(qtl_dat)\n",
" context = \"${_meta[0]}\" %>% str_split(.,\"@\",simplify = T) %>% .[,1] \n",
" gene_region = pecotmr:::get_nested_element(qtl_dat[[1]], c(context,\"region_info\", \"grange\"))\n",
" \n",
" # Step 1: find relevant GWAS regions that overlap each the xQTL region of interest\n",
" gwas_finemapped_datas <- c(${ \",\".join(['\"%s\"' % x for x in _meta[1:]])}) \n",
" gwas_finemapped_datas <- gwas_finemapped_datas[grep('rds$',gwas_finemapped_datas)]\n",
" \n",
" gwas_regions <- gwas_finemapped_datas %>% basename %>% str_extract(., \"chr\\\\d+_\\\\d+_\\\\d+\")\n",
" overlap_index <- NULL\n",
" for (i in 1:length(gwas_regions)) {\n",
" region <- gwas_regions[i]\n",
" split_region <- unlist(strsplit(region, \"_\"))\n",
" block_chrom <- as.numeric(split_region[1] %>% gsub(\"chr\",\"\",.))\n",
" block_start <- as.numeric(split_region[2] %>% strsplit(., \"-\") %>% unlist %>% .[1])\n",
" block_end <- as.numeric(split_region[2] %>% strsplit(., \"-\") %>% unlist %>% .[2])\n",
" if (gene_region$chrom == block_chrom && (gene_region$start <= block_end | gene_region$end >= block_start)) {\n",
" overlap_index <- c(overlap_index, i)\n",
" }\n",
" }\n",
" \n",
" if (!is.null(overlap_index)) {\n",
" gwas_finemapped_data <- gwas_finemapped_datas[overlap_index]\n",
" \n",
" # Step 2: load enrichment analysis results\n",
" # Extract values for p1, p2, and p12\n",
" if(${\"FALSE\" if skip_enrich else \"TRUE\"}){\n",
" enrich_file = paste0('${cwd:a}','/', '${name}', '.' ,context,'.enrichment.rds')\n",
" p1 <- readRDS(enrich_file)[[1]]$`Alternative (coloc) p1`\n",
" p2 <- readRDS(enrich_file)[[1]]$`Alternative (coloc) p2`\n",
" p12 <- readRDS(enrich_file)[[1]]$`Alternative (coloc) p12`\n",
" \n",
" message(\"Priors are P1:\", p1, \"; p2: \", p2, \"; p12: \", p12)\n",
" } else {\n",
" p1 = 1e-4\n",
" p2 = 1e-4\n",
" p12 = 5e-6\n",
" message(\"Priors are using default, P1: 1e-4, P2: 1e-4, P12: 5e-6\" )\n",
" }\n",
" # Step 3: Apply colocalization analysis between each condition and GWAS\n",
" coloc_res[[gene]] <- coloc_wrapper(xqtl_file = xqtl_finemapped_data, gwas_files = gwas_finemapped_data, \n",
" xqtl_finemapping_obj = c(context,${\",\".join(['\"%s\"' % x for x in xqtl_finemapping_obj]) if len(xqtl_finemapping_obj) != 0 else \"NULL\"}), \n",
" xqtl_varname_obj = c(context,${\",\".join(['\"%s\"' % x for x in xqtl_varname_obj]) if len(xqtl_varname_obj) != 0 else \"NULL\"}), \n",
" gwas_finemapping_obj = c(${\",\".join(['\"%s\"' % x for x in gwas_finemapping_obj]) if len(gwas_finemapping_obj) != 0 else \"NULL\"}), \n",
" gwas_varname_obj = c(${\",\".join(['\"%s\"' % x for x in gwas_varname_obj]) if len(gwas_varname_obj) != 0 else \"NULL\"}),\n",
" xqtl_region_obj = c(context,${\",\".join(['\"%s\"' % x for x in xqtl_region_obj]) if len(xqtl_region_obj) != 0 else \"NULL\"}), \n",
" gwas_region_obj = c(${\",\".join(['\"%s\"' % x for x in gwas_region_obj]) if len(gwas_region_obj) != 0 else \"NULL\"}),\n",
" p1 = p1, p2 = p2, p12 = p12, filter_lbf_cs = ${\"FALSE\" if skip_enrich else \"TRUE\"})\n",
"\n",
" if (${\"TRUE\" if ld_meta_file_path.is_file() else \"FALSE\"}) {\n",
" if(length(coloc_res[[gene]]) > 2) {\n",
" coloc_res[[gene]] <- coloc_post_processor(coloc_res[[gene]], LD_meta_file_path = ${ld_meta_file_path:r}, analysis_region = coloc_res[[gene]]$analysis_region)\n",
" if(!is.null(coloc_res[[gene]]$sets$cs)) writeLines(coloc_res[[gene]]$sets$cs %>% unlist, gsub(\"rds$\",\"coloc_res\",\"${_output}\"))\n",
" } else { message(\"No coloc results were generated.\") }\n",
" }\n",
" \n",
" } else {\n",
" message(\"No overlap was found between GWAS blocks and QTL analysis region.\")\n",
" coloc_res <- \"No overlap was found between GWAS blocks and QTL analysis region.\"\n",
" }\n",
" }\n",
" saveRDS(coloc_res, ${_output:r})"
"bash: expand = '${ }', stdout = f\"{_output:n}.stdout\", stderr = f\"{_output:n}.stderr\", container = container\n",
" set -e\n",
" # context (before '@' in _meta[0]) selects the per-region enrichment file the\n",
" # enrichment step wrote: {cwd}/{name}.{context}.enrichment.rds (legacy naming).\n",
" # Convert this region's xQTL + GWAS sides, then run colocPipeline (enloc\n",
" # variant when --enrichment is supplied).\n",
" Rscript code/script/pecotmr_integration/legacy_enloc_finemap_convert.R \\\n",
" --mode qtl \\\n",
" --rds-files ${\",\".join([str(x) for x in _input])} \\\n",
" --finemapping-obj '${\" \".join(xqtl_finemapping_obj)}' \\\n",
" --varname-obj '${\" \".join(xqtl_varname_obj)}' \\\n",
" --study ${name} \\\n",
" --output ${_output:n}.qtl.rds\n",
"\n",
" Rscript code/script/pecotmr_integration/legacy_enloc_finemap_convert.R \\\n",
" --mode gwas \\\n",
" --rds-files ${\",\".join([str(x) for x in _meta[1:] if str(x).endswith(\".rds\")])} \\\n",
" --finemapping-obj '${\" \".join(gwas_finemapping_obj)}' \\\n",
" --varname-obj '${\" \".join(gwas_varname_obj)}' \\\n",
" --output ${_output:n}.gwas.rds\n",
"\n",
" Rscript code/script/pecotmr_integration/coloc.R \\\n",
" --qtl-fine-mapping ${_output:n}.qtl.rds \\\n",
" --gwas-input ${_output:n}.gwas.rds \\\n",
" ${('--enrichment '+cwd.absolute().joinpath(name + \".\" + _meta[0].split(\"@\")[0] + \".enrichment.rds\").as_posix()) if not skip_enrich else ''} \\\n",
" --output ${_output}\n"
]
},
{
Expand Down Expand Up @@ -722,4 +687,4 @@
},
"nbformat": 4,
"nbformat_minor": 5
}
}
Loading
Loading