xQTL-GWAS pairwise enrichment and colocalization

sos run pipeline/SuSiE_enloc.ipynb xqtl_gwas_enrichment \
    --gwas-meta-data tests/fixtures/susie_enloc/protocol_example.enloc.gwas_meta.tsv \
    --xqtl-meta-data tests/fixtures/susie_enloc/protocol_example.enloc.xqtl_meta.tsv \
    --xqtl-finemapping-obj preset_variants_result susie_result_trimmed \
    --xqtl-varname-obj preset_variants_result variant_names \
    --gwas-finemapping-obj AD_Bellenguez_2022 RSS_QC_RAISS_imputed susie_result_trimmed \
    --gwas-varname-obj  AD_Bellenguez_2022 RSS_QC_RAISS_imputed variant_names \
    --xqtl-region-obj   region_info grange \
    --qtl-path output/xqtl_gwas_enrichment/susie_enloc_data \
    --gwas-path output/xqtl_gwas_enrichment/susie_enloc_data \
    --context-meta tests/fixtures/susie_enloc/protocol_example.enloc.context_meta.tsv \
    --cwd output/xqtl_gwas_enrichment

sos run pipeline/SuSiE_enloc.ipynb susie_coloc \
    --gwas-meta-data tests/fixtures/susie_enloc/protocol_example.enloc.gwas_meta.tsv \
    --xqtl-meta-data tests/fixtures/susie_enloc/protocol_example.enloc.xqtl_meta.tsv \
    --xqtl-finemapping-obj preset_variants_result susie_result_trimmed \
    --xqtl-varname-obj preset_variants_result variant_names \
    --gwas-finemapping-obj AD_Bellenguez_2022 RSS_QC_RAISS_imputed susie_result_trimmed \
    --gwas-varname-obj  AD_Bellenguez_2022 RSS_QC_RAISS_imputed variant_names \
    --xqtl-region-obj   region_info grange \
    --qtl-path output/xqtl_gwas_enrichment/susie_enloc_data \
    --gwas-path output/xqtl_gwas_enrichment/susie_enloc_data \
    --context-meta tests/fixtures/susie_enloc/protocol_example.enloc.context_meta.tsv \
    --ld-meta-file-path tests/fixtures/ld_reference/ld_meta_file.tsv \
    --skip-enrich \
    --cwd output/susie_coloc

sos run pipeline/SuSiE_enloc.ipynb -h

[global]
# Workdir
parameter: cwd = path("output")
parameter: modular_script_dir = path('code/script')
# A list of file paths for fine-mapped GWAS results. 
parameter: gwas_meta_data = path
# A list of file paths for fine-mapped xQTL results. 
parameter: xqtl_meta_data = path
# Optional: if a region list is provide the enrichment analysis will be focused on provided region. 
# The LAST column of this list will contain the ID of regions to focus on
parameter: region_list = path()
# Optional: if a region name is provided 
# the analysis would be focused on the union of provides region list and region names
parameter: region_name = []
# It is required to input the name of the analysis
parameter: name = f"{xqtl_meta_data:bnn}.{gwas_meta_data:bnn}"
parameter: container = ""
# For cluster jobs, number commands to run per job
parameter: job_size = 200
# Wall clock time expected
parameter: walltime = "5m"
# Memory expected: quite large for enrichment analysis but small for xQTL colocalization
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 1
#Optional table name in xQTL RDS files to get fine mapping results (eg 'susie_result_trimmed ').
parameter: xqtl_finemapping_obj = []
#Optional table name in xQTL RDS files to get variable names (eg ' variant_names ').
parameter: xqtl_varname_obj = []
#Optional table name in GWAS RDS files to get fine mapping results (eg 'AD_Bellenguez_2022 single_effect_regression susie_result_trimmed ').
parameter: gwas_finemapping_obj = []
#Optional table name in GWAS RDS files to get variable names(eg 'AD_Bellenguez_2022 single_effect_regression variant_names ').
parameter: gwas_varname_obj = []
#Optional table name in xQTL RDS files to get region info (eg 'susie_result_trimmed region_info grange ').
parameter: xqtl_region_obj = []
#Optional table name in GWAS RDS files to get region info (eg 'AD_Bellenguez_2022 single_effect_regression region_info grange ').
parameter: gwas_region_obj = []
#Directory path for GWAS orignal finemapping results 
parameter: gwas_path = ''
#Directory path for xQTL orignal finemapping results 
parameter: qtl_path = ''
# a meta file showing the context and corresponding analysis_name
parameter: context_meta = path()
# Optional: if a region list is provide the analysis will be focused on provided region. 
# The LAST column of this list will contain the ID of regions to focus on
# Otherwise, all regions with both genotype and phenotype files will be analyzed
parameter: region_list = path()
# Optional: if a region name is provided 
# the analysis would be focused on the union of provides region list and region names
parameter: region_name = []
# use conditions_top_loci_minp column or not, which can reduce a lot computing resource by pick top 1 isoform
parameter: minp = False

[get_analysis_regions: shared = {'enrichment_rows': "list(__import__('csv').DictReader(open(enrichment_manifest), delimiter=chr(9)))"}]
# Resolve the per-condition xQTL x GWAS enrichment units into a manifest TSV via
# enloc_manifest.R -- the coordinate-overlap pairing + per-condition study
# routing that the legacy get_analysis_regions did in-notebook with pandas. The
# parsed rows are shared as enrichment_rows, so [xqtl_gwas_enrichment]
# auto-triggers this step via depends: sos_variable -- the notebook is invoked
# exactly as before (e.g. sos run SuSiE_enloc.ipynb xqtl_gwas_enrichment).
input: None
enrichment_manifest = f"{cwd:a}/get_analysis_regions/{name}.enrichment_manifest.tsv"
bash: expand = "${ }", container = container
    mkdir -p ${cwd:a}/get_analysis_regions
    Rscript ${modular_script_dir}/pecotmr_integration/enloc_manifest.R --mode enrichment \
        --xqtl-meta ${xqtl_meta_data} --gwas-meta ${gwas_meta_data} \
        ${("--context-meta " + str(context_meta)) if context_meta.is_file() else ""} \
        ${("--qtl-path " + str(qtl_path)) if qtl_path else ""} \
        ${("--gwas-path " + str(gwas_path)) if gwas_path else ""} \
        ${("--region-list " + str(region_list)) if region_list.is_file() else ""} \
        ${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
        ${("--gwas-finemapping-obj " + " ".join([str(x) for x in gwas_finemapping_obj])) if len(gwas_finemapping_obj) > 0 else ""} \
        ${"--minp" if minp else ""} \
        --output ${enrichment_manifest}

[get_overlapped_analysis_regions: shared = {'coloc_rows': "list(__import__('csv').DictReader(open(coloc_manifest), delimiter=chr(9)))"}]
# Resolve the per-(condition,region) xQTL x GWAS coloc units into a manifest TSV
# via enloc_manifest.R --mode coloc -- the block_top_loci variant-overlap pairing
# + per-condition study routing that the legacy get_overlapped_analysis_regions
# did in-notebook with pandas. Shared as coloc_rows, so [susie_coloc]
# auto-triggers this step (sos run SuSiE_enloc.ipynb susie_coloc) unchanged.
input: None
coloc_manifest = f"{cwd:a}/get_overlapped_analysis_regions/{name}.coloc_manifest.tsv"
bash: expand = "${ }", container = container
    mkdir -p ${cwd:a}/get_overlapped_analysis_regions
    Rscript ${modular_script_dir}/pecotmr_integration/enloc_manifest.R --mode coloc \
        --xqtl-meta ${xqtl_meta_data} --gwas-meta ${gwas_meta_data} \
        ${("--context-meta " + str(context_meta)) if context_meta.is_file() else ""} \
        ${("--qtl-path " + str(qtl_path)) if qtl_path else ""} \
        ${("--gwas-path " + str(gwas_path)) if gwas_path else ""} \
        ${("--region-list " + str(region_list)) if region_list.is_file() else ""} \
        ${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
        ${("--gwas-finemapping-obj " + " ".join([str(x) for x in gwas_finemapping_obj])) if len(gwas_finemapping_obj) > 0 else ""} \
        ${"--minp" if minp else ""} \
        --output ${coloc_manifest}

[xqtl_gwas_enrichment]
depends: sos_variable("enrichment_rows")
# Reproducibility: integer RNG seed for the enrichment multiple-imputation
# sampler, which draws in C++ and so is unreachable from R's set.seed; -1 = unset
# (nondeterministic).
parameter: seed = -1
seed_arg = f"--seed {seed}" if seed >= 0 else ""
jobs = enrichment_rows
stop_if(len(jobs) == 0, f'No files left for analysis')
_qtl_groups = [j['qtl_files'].split(',') for j in jobs]
input: [f for grp in _qtl_groups for f in grp], group_by = lambda x: _qtl_groups, group_with = "jobs"
output: f'{cwd:a}/{name}.{_jobs["unit_id"]}.enrichment.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = '${ }', stdout = f"{_output:n}.stdout", stderr = f"{_output:n}.stderr", container = container
    set -e
    # The xQTL and GWAS sides are already S4 FineMappingResult collections
    # (fine_mapping output); qtl_enrichment.R combines its multi-file inputs
    # itself, so the legacy legacy_enloc_finemap_convert.R step is gone.
    Rscript ${modular_script_dir}/pecotmr_integration/qtl_enrichment.R \
        --qtl-fine-mapping ${" ".join([str(x) for x in _input])} \
        --gwas-fine-mapping ${_jobs["gwas_files"].replace(",", " ")} \
        --ncore ${numThreads} \
        ${seed_arg} \
        --output ${_output}

[susie_coloc]
depends: sos_variable("coloc_rows")
# depends: sos_step("xqtl_gwas_enrichment") #changed
stop_if(len(coloc_rows) == 0, f'No files left for analysis')
# skip enrichment or not
parameter: skip_enrich=False
# filter lbf by cs as default in susie coloc. default is False means filtering by V > prior_tol
parameter: filter_lbf_cs = False
# ld reference path for postprocessing 
parameter: ld_meta_file_path = path()
jobs = coloc_rows
_qtl_groups = [j['qtl_files'].split(',') for j in jobs]
input: [f for grp in _qtl_groups for f in grp], group_by = lambda x: _qtl_groups, group_with = "jobs"
output: f'{cwd:a}/{step_name}/{name}.{_jobs["unit_id"]}.coloc.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = '${ }', stdout = f"{_output:n}.stdout", stderr = f"{_output:n}.stderr", container = container
    set -e
    # xQTL + GWAS sides are already S4 FineMappingResult collections; coloc.R
    # combines its multi-file inputs itself (no legacy conversion). The enrichment
    # file (context before '@' in unit_id) tunes colocPipeline's prior.
    Rscript ${modular_script_dir}/pecotmr_integration/coloc.R \
        --qtl-fine-mapping ${" ".join([str(x) for x in _input])} \
        --gwas-input ${_jobs["gwas_files"].replace(",", " ")} \
        ${('--filter-lbf-cs' if filter_lbf_cs else '')} \
        ${('--enrichment '+cwd.absolute().joinpath(name + "." + _jobs["unit_id"].split("@")[0] + ".enrichment.rds").as_posix()) if not skip_enrich else ''} \
        --output ${_output}