Fine-mapping result post-processing

Turns per-region fine-mapping objects into flat TSV tables, diagnostic plots, VCFs and QTL-GWAS overlap summaries.

Overview

Fine-mapping leaves you with per-region R objects, which are the wrong shape for almost everything that comes next - inspecting results, sharing them, or intersecting them with a GWAS. This module reads those objects and re-expresses them as text: susie_to_tsv writes a per-variant table (PIP, credible-set membership, posterior mean and sd), a wide per-variant by per-effect log-Bayes-factor table, and a per-credible-set effect table. The remaining workflows build on those tables - collapsing them across regions, drawing PIP landscape and UpSet plots, annotating variants, exporting VCFs, and overlapping QTL credible sets with GWAS signals.

The export step consumes a FineMappingResult object rather than a raw susieR fit, so the three tables are views onto that object: variant.tsv is its topLoci view, lbf.tsv its lbf view and effect.tsv its cs_summary view. Because the object carries its own region and trait, --region-list is accepted for CLI compatibility but is no longer used.

When to run it. After fine-mapping, on the outputs of susie_twas or fsusie, and before any downstream reporting or GWAS integration. Most of the workflows beyond susie_to_tsv additionally need GWAS summary statistics, plot-recipe files or a singularity container.

Input

• --rds_path input/mnm_postprocessing/protocol_example.ENSG00000283047.fine_mapping.rds (one or more FineMappingResult RDS files written by mnm_regression; each carries its own region and trait, and the step runs once per file)

QtlFineMappingResult: S4 DataFrame, 2 rows x 7 columns
$ study: character [1:2] test_study test_study
$ context: character [1:2] context1 context2
$ trait: character [1:2] ENSG00000283047 ENSG00000283047
$ method: character [1:2] susie susie
$ entry: S4 object  <SimpleList>
 ..$ listData: List of 2
 .. ..$ : S4 object  <FineMappingEntry>
 .. .. ..$ variantIds: ...
 .. .. ..$ susieFit: ...
 .. .. ..$ topLoci: ...
 .. .. ..$ cvResult: ...
 .. ..$ : S4 object  <FineMappingEntry>
 .. .. ..$ variantIds: ...
 .. .. ..$ susieFit: ...
 .. .. ..$ topLoci: ...
 .. .. ..$ cvResult: ...
 ..$ elementType: character [1:1] ANY
 ..$ elementMetadata: name [1:1] 
NULL

 ..$ metadata: List of 0
... (truncated)

• --study toy (tag used as name in output filenames and in intermediate RDS names. Required.)

• --cwd output_postproc (working directory all outputs are written under. Defaults to output.)

• --container (an empty string runs the step natively in the project environment instead of through a singularity image)

• --region-list (accepted for CLI compatibility and ignored. The FineMappingResult carries its own region and trait, so no external region list is needed.)

The same directory also holds protocol_example.fsusie.fine_mapping.rds and protocol_example.mvsusie.fine_mapping.rds, for the fsusie_to_tsv and mv_susie_2 workflows.

The export_top_loci workflow consolidates results across regions and takes a different set of inputs, none of which the example data ships:

• --file_path (directory holding the per-region fine-mapping RDS files)

• --prefix, --suffix (the filename parts before and after the region id, so each region resolves to <file_path>/<prefix>.<region>.<suffix> - for example ROSMAP_AC_eQTL.ENSG00000012779.univariate_bvsr.rds)

• --qtl_type cis (one of cis, trans or trans_<program>; appears in the combined output filename)

• --region_file (BED file listing the regions to consolidate)

The remaining workflows take these additional flags:

• --tsv_path output_postproc/*.lbf.tsv (per-region LBF tables to concatenate; susie_tsv_collapse)

• --annot_tibble <annotation_tibble.tsv> (annotation tibble used to label variants; susie_pip_landscape_plot and fsusie_extract_effect)

• --trait_to_select 1 (how many traits the credible-set UpSet diagram shows; susie_upsetR_cs_plot)

• --SNP_list <variant_list.txt> (variants to annotate; tmp_annotatation_of_snps)

• --annotation_rds <annotation.rds> (annotation reference the variants are looked up against; tmp_annotatation_of_snps)

• --geno_ref <genotype_reference.bim.gz> (genotype reference used to resolve variant identifiers; export_top_loci)

• --min_corr 0.8 (minimum correlation for a variant to join a top-loci block; export_top_loci)

• --condition_meta <gwas_meta.tsv> (GWAS meta file naming the conditions to overlap; overlap_qtl_gwas)

Output

Per-region TSV export (susie_to_tsv, sQTL_susie_to_tsv, fsusie_to_tsv), one set per input SuSiE object:

• <region>.variant.tsv – one row per variant: PIP, credible-set order and id, posterior mean and sd.

• <region>.lbf.tsv – per-variant log Bayes factors for each single-effect.

• <region>.effect.tsv – per-effect summary (susie_to_tsv and sQTL_susie_to_tsv only).

• <name>.all_variants.tsv – variants pooled across regions (fsusie_to_tsv only).

Collapsed tables (susie_tsv_collapse):

• <prefix>.chr<chromosome>.unisusie_rss.lbf.tsv – per-region lbf tables concatenated by chromosome.

Diagnostic plots:

• <region>.pip_landscape_plot.pdf (susie_pip_landscape_plot)

• <name>.UpSetR.pdf and <name>.UpSetR_cs.pdf (susie_upsetR_plot, susie_upsetR_cs_plot)

Variant annotation (tmp_annotatation_of_snps):

• <name>.annotated.rds and <name>.annotated_rev.rds

fSuSiE effect extraction:

• <region>.estimated_effect.tsv (fsusie_extract_effect)

• <region>.affected_region.tsv (fsusie_affected_region)

VCF export (mv_susie_2):

• <region>.<context>.vcf.bgz

cis and GWAS consolidation:

• <name>_cache/<name>_<meta_info>.tsv, collapsed into <name>.cis_results_db.tsv (or <name>.block_results_db.tsv for GWAS input) – cis_results_export

• <name>.union_export.tsv.gz – gwas_results_export

• summary/<region>.top_loci.bed.gz, collapsed into summary/<name>.<qtl_type>.top_loci.bed.gz – export_top_loci

QTL and GWAS overlap (overlap_qtl_gwas):

• gwas_qtl/cache/<name>_gwas_batch_meta_<meta_info>.tsv

• <name>.overlapped.gwas.tsv

No example output ships with the repository: output_postproc/ is empty and the pipeline cannot be re-run here, so no structure dumps are shown. The filenames above come from each step output: statement.

Minimal Working Example

Export fine-mapping results to per-variant, LBF and effect tables

susie_to_tsv reads each FineMappingResult object and writes its three text views.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb susie_to_tsv \
    --cwd output_postproc --study toy \
    --rds_path input/mnm_postprocessing/protocol_example.ENSG00000283047.fine_mapping.rds \
    --container '' -j1

fSuSiE and sQTL exports

fsusie_to_tsv and sQTL_susie_to_tsv perform the same export for fSuSiE and sQTL objects.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb fsusie_to_tsv \
    --cwd output_postproc --study <study_name> \
    --rds_path input/mnm_postprocessing/protocol_example.fsusie.fine_mapping.rds \
    --container '' -j1

Collapse per-region tables

susie_tsv_collapse concatenates the per-region LBF tables into one table per chromosome.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb susie_tsv_collapse \
    --cwd output_postproc --study <study_name> \
    --tsv_path output_postproc/*.lbf.tsv \
    --container '' -j1

Diagnostic plots

susie_pip_landscape_plot draws the PIP landscape; susie_upsetR_plot and susie_upsetR_cs_plot draw UpSet diagrams of credible-set overlap. The landscape plot needs an annotation tibble; --trait_to_select picks how many traits the credible-set diagram shows.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb susie_pip_landscape_plot \
    --cwd output_postproc --study <study_name> \
    --rds_path <fine_mapping_dir>/*.fine_mapping.rds \
    --annot_tibble <annotation_tibble.tsv> \
    --container '' -j1

sos run pipeline/mnm_postprocessing.ipynb susie_upsetR_cs_plot \
    --cwd output_postproc --study <study_name> \
    --rds_path <fine_mapping_dir>/*.fine_mapping.rds \
    --trait_to_select 1 \
    --container '' -j1

Variant annotation

tmp_annotatation_of_snps annotates a variant list against a reference RDS, then writes the reverse mapping.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb tmp_annotatation_of_snps \
    --cwd output_postproc --study <study_name> \
    --SNP_list <variant_list.txt> \
    --annotation_rds <annotation.rds> \
    --container '' -j1

fSuSiE effect extraction

fsusie_extract_effect writes the estimated effect curve per region; fsusie_affected_region writes the intervals a signal spans.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb fsusie_extract_effect \
    --cwd output_postproc --study <study_name> \
    --rds_path input/mnm_postprocessing/protocol_example.fsusie.fine_mapping.rds \
    --annot_tibble <annotation_tibble.tsv> \
    --container '' -j1

VCF export

mv_susie_2 writes multivariate results as bgzipped VCF, one file per region and context.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb mv_susie \
    --cwd output_postproc --study <study_name> \
    --rds_path input/mnm_postprocessing/protocol_example.mvsusie.fine_mapping.rds \
    --container '' -j1

cis and GWAS consolidation

cis_results_export collects per-region results into a results database, gwas_results_export does the same for GWAS blocks, combine_export_meta merges the per-chunk caches, and export_top_loci writes the combined top-loci BED. All four resolve each region as <file_path>/<prefix>.<region>.<suffix> and need a GWAS meta file.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb export_top_loci \
    --cwd output_postproc --study <study_name> \
    --region_file <regions.bed> \
    --file_path <fine_mapping_dir> \
    --prefix <fm_rds_prefix> --suffix <fm_rds_suffix> \
    --qtl_type cis --min_corr 0.8 \
    --geno_ref <genotype_reference.bim.gz> \
    --container '' -j1

QTL and GWAS overlap

overlap_qtl_gwas intersects QTL credible sets with GWAS results, caching a batch meta file per chunk before writing the overlap table.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_postprocessing.ipynb overlap_qtl_gwas \
    --cwd output_postproc --study <study_name> \
    --region_file <regions.bed> \
    --file_path <fine_mapping_dir> \
    --prefix <fm_rds_prefix> --suffix <fm_rds_suffix> \
    --condition_meta <gwas_meta.tsv> \
    --container '' -j1

The meta file consumed by the consolidation workflows is produced by:

get_condition <- function(conditions, Author, qtl_type){
    strings <- c()
    for(condition in conditions){
        string =  paste(condition, Author, qtl_type, sep = "_")  
        string = paste(unique(unlist(strsplit(string, "_"))), collapse = "_")
        strings <- c(strings, string)
    }
    return(strings)
}

raw_name<- c("Mic","Ast","Oli","OPC","Exc","Inh","DLPFC","PCC","AC")
raw_name_kellis<- c("Mic_Kellis","Ast_Kellis","Oli_Kellis","OPC_Kellis","Exc_Kellis","Inh_Kellis","Ast.10","Mic.12","Mic.13")

dejager_name <- get_condition(raw_name, "De_Jager","eQTL")
kellis_name <- get_condition(raw_name_kellis, "Kellis","eQTL")
eQTL_meta <- data.frame(raw_name = c(raw_name, raw_name_kellis), new_name = c(dejager_name, kellis_name))
write_delim(eQTL_meta, "output/mnm_postprocessing/eQTL_meta.tsv")

Command Interface

sos run pipeline/mnm_postprocessing.ipynb -h

usage: sos run pipeline/mnm_postprocessing.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  susie_to_tsv
  sQTL_susie_to_tsv
  fsusie_to_tsv
  susie_tsv_collapse
  susie_pip_landscape_plot
  susie_upsetR_plot
  susie_upsetR_cs_plot
  tmp_annotatation_of_snps
  fsusie_extract_effect
  fsusie_affected_region
  mv_susie
  cis_results_export
  gwas_results_export
  export_top_loci
  combine_export_meta
  overlap_qtl_gwas

Global Workflow Options:
  --cwd output (as path)
                        A region list file documenting the chr_pos_ref_alt of a
                        susie_object
  --study VAL (as str, required)
  --container ''
                        Path to work directory where output locates Containers
                        that contains the necessary packages
  --modular-script-dir code/script (as path)
  --entrypoint  ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""

  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 96h
                        Wall clock time expected
  --mem 6G
                        Memory expected
  --numThreads 2 (as int)
                        Number of threads
  --windows 1000000 (as int)

Sections
  susie_to_tsv_1:
    Workflow Options:
      --region-list . (as path)
                        Migrated: consumes a FineMappingResult RDS
                        (fine_mapping.R output) instead of a raw susieR object.
                        The three legacy outputs map to fine_mapping_export.R
                        views: variant.tsv = topLoci   (per-variant: pip,
                        posterior mean/sd, logBF, cs, region) lbf.tsv     = lbf
                        (wide per-variant x per-effect log Bayes factor matrix)
                        effect.tsv  = cs_summary (per credible set: size,
                        purity, V, logBF, lead variant) The FineMappingResult
                        carries its own region/trait, so `region_list` is
                        accepted for CLI compatibility but is no longer used.
      --rds-path  paths

  sQTL_susie_to_tsv_1:
    Workflow Options:
      --region-list . (as path)
                        Migrated: identical to susie_to_tsv but preserves the
                        sQTL / leafcutter2 filename sanitization ("*" -> "N";
                        ":"/"+" scrubbed from the task tag), since intron ids
                        contain shell-hostile characters. Consumes a
                        FineMappingResult RDS and writes the topLoci / lbf /
                        cs_summary views.
      --rds-path  paths

  fsusie_to_tsv_1:
    Workflow Options:
      --region-list . (as path)
                        Migrated: consumes an fSuSiE FineMappingResult RDS.
                        variant.tsv = topLoci (per-variant: pip, cs,
                        cs_95_purity) lbf.tsv     = lbf     (fSuSiE per-variant
                        x per-effect lBF matrix) (Fixes the legacy bug where
                        lbf.tsv was written the variant table, not the lbf
                        matrix.) is_dummy = cs_95_purity < 0.5; the functional
                        effect curves + peak positions are exported by
                        fsusie_extract_effect / fsusie_affected_region
                        (credible_band / affected_regions views), so they are
                        not duplicated here.
      --rds-path  paths

  *_to_tsv_2:
    Workflow Options:
      --name  f'{_input[0]:b}'.split(".")[0]

                        Reduce: concatenate the per-region variant tables (the
                        *variant.tsv outputs of the parallel _1 map step) into
                        one file, in R via data.table::rbindlist (fill=TRUE
                        aligns columns). _input carries all of _1's outputs, so
                        keep only the variant tables (not the lbf / effect
                        ones).
  susie_tsv_collapse:
    Workflow Options:
      --tsv-path  paths

                        Concatenate the per-region lbf tables (named
                        *.chr<N>_<start>_<end>...lbf.tsv) by chromosome.
                        Chromosome split is plain python (no pandas); the concat
                        is R data.table::rbindlist (headers handled by fread;
                        fill=TRUE aligns columns).
  susie_pip_landscape_plot:
    Workflow Options:
      --rds-path  paths

                        Per-region PIP landscape from a FineMappingResult
                        (fine_mapping output). Modernized from the legacy multi-
                        TSV plot: pecotmr does no plotting, so the ggplot lives
                        in fine_mapping_pip_plot.R. (The regulatory/gene
                        annotation overlay is a deferred enhancement;
                        annot_tibble is kept for CLI compat.)
  susie_upsetR_plot:
    Workflow Options:
      --rds-path  paths

                        UpSet plot of credible-set variant overlap across the
                        FineMappingResult(s)' (study, context, trait, method)
                        tuples. pecotmr does no plotting; the UpSet lives in
                        fine_mapping_upset.R.
  susie_upsetR_cs_plot:
    Workflow Options:
      --rds-path  paths

                        Credible-set overlap UpSet plot across the
                        FineMappingResult(s), via fine_mapping_upset.R (pecotmr
                        does no plotting). trait_to_select is retained for CLI
                        compatibility; the single-trait "CS shared with other
                        phenotypes" sub-view is a deferred wrapper option.
      --trait-to-select 1 (as int)
  tmp_annotatation_of_snps_1:
    Workflow Options:
      --SNP-list VAL (as path, required)
      --annotation-rds VAL (as path, required)
  tmp_annotatation_of_snps_2:
    Workflow Options:
      --SNP-list VAL (as path, required)
      --annotation-rds VAL (as path, required)
  fsusie_extract_effect:
    Workflow Options:
      --rds-path  paths

                        Fitted effect curve + credible band per fSuSiE credible
                        set, via fsusieCredibleBand (fine_mapping_export.R
                        `credible_band` view). pecotmr does the wavelet credible
                        band; it needs an UNtrimmed fSuSiE fit (fine_mapping
                        trim=FALSE). The effect-curve PDF is a deferred plot
                        wrapper (pecotmr does no plotting); annot_tibble is kept
                        for CLI compatibility.
  fsusie_affected_region:
    Workflow Options:
      --rds-path  paths

                        GRanges of intervals where the fSuSiE credible band
                        excludes 0 (with a direction mcol), via
                        fsusieAffectedRegions (fine_mapping_export.R
                        `affected_regions` view). Needs an UNtrimmed fSuSiE fit
                        (fine_mapping trim=FALSE). The region PDF is a deferred
                        plot wrapper.
  mv_susie_2:
    Workflow Options:
      --rds-path  paths

                        Write per-context fine-mapping VCFs (ES / SE / LP / AF +
                        PIP / CS / within_cs_pip in the sample column) from a
                        FineMappingResult via fine_mapping_vcf.R ->
                        pecotmr::writeSumstatsVcf, replacing the inline
                        create_vcf. writeSumstatsVcf splits a multi-context
                        collection into one .vcf.bgz per context; the runtime
                        set is tracked with output: dynamic() (so {cwd} must be
                        created up front, since a dynamic output -- unlike a
                        concrete target -- does not pre-create its dir).
  cis_results_export_1, gwas_results_export_1, export_top_loci_1: Exporting cis
                        susie_twas results
    Workflow Options:
      --per-chunk 200 (as int)
                        per chunk we process at most 200 datasets
      --region-file . (as path)
                        Region list should have last column being region name.
      --file-path ''
                        the path stored original output files
      --prefix  (as list)
                        assuming orinal files are named as prefix.id.suffix
                        prefix of the original files (before id) to identify
                        file.
      --suffix VAL (as str, required)
                        suffix of the original files (after id) to identify
                        file.
      --condition-meta . (as path)
                        if need to rename context, please load condition meta
                        file
      --pip-thres 0.05 (as float)
                        if the pip all variants in a cs < this threshold, then
                        remove this cs
      --cs-size 3 (as int)
                        only keep top cs_size variants in one cs
      --exported-file . (as path)
                        provide exported meta to filter the exported genes
      --region-list . (as 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
      --region-name  (as list)
                        Optional: if a region name is provided the analysis
                        would be focused on the union of provides region list
                        and region names
      --geno-ref . (as path)
                        the (aligned) geno bim file to check allele flipping.
                        Retained for CLI compatibility; allele alignment is now
                        done upstream by fine_mapping.R (matchVariants), so it
                        is unused here.
      --context-meta . (as path)
                        context meta file to map renamed context backs
      --min-corr  (as list)
                        CS purity (min.abs.corr) cutoff(s); the first value is
                        used for conditions_top_loci.
      --[no-]gwas (default to False)
                        set this parameter as `True` when exporting gwas data.
      --[no-]fsusie (default to False)
                        set this parameter as `True` when exporting fsusie data.
      --[no-]metaQTL (default to False)
                        set this parameter as `True` when exporting metaQTL
                        data.
      --[no-]mnm (default to False)
                        set this parameter as `True` when exporting mnm data.
  cis_results_export_2, gwas_results_export_2:
    Workflow Options:
      --[no-]step1-only (default to False)
                        set it as True if you only want run first step (running
                        parallel) and combine them manually after all finished
      --exported-file . (as path)
                        provide exported meta to filter the exported genes
      --[no-]gwas (default to False)
                        optional: qtl or gwas, there is slightly different in
                        qtl and gwas rds file
  cis_results_export_3:
    Workflow Options:
      --[no-]step1-only (default to False)
                        set it as True if you only want run first step (running
                        parallel) and combine them manually after all finished
  gwas_results_export_3: get union of step1 1200 blocks costed ~5mins with in
                        one for loop
    Workflow Options:
      --[no-]step1-only (default to False)
  combine_export_meta:  simply combine seperate meta files, not works for having
                        an exsiting exported file for now.
    Workflow Options:
      --cache-path VAL (as path, required)
      --output-file VAL (as str, required)
      --[no-]remove-cache (default to False)
  export_top_loci_2:
    Workflow Options:
      --export-path  cwd

                        Export the per-region top-loci BED from each
                        cis_results_db (an FMR in the modern pipeline) via the
                        pecotmr-backed fine_mapping_top_loci_bed.R wrapper.
                        getTopLoci does the PIP/CS + independent purity
                        filtering; lfsr / conditional_effect are per-(variant,
                        context) numeric columns, so the legacy filter_lfsr
                        semicolon split/rejoin is gone. `export_suffix`
                        deliberately differs from step 1's `--suffix` (input
                        rds) to avoid CLI collision.
      --export-prefix ''
      --export-suffix 'cis_results_db.rds'
      --fsusie-prefix ''
      --[no-]preset-top-loci (default to False)
      --signal-cutoff 0.025 (as float)
                        PIP cutoff for getTopLoci; optional independent CS
                        purity (min.abs.corr) cutoff (<0 => none).
      --min-purity -1.0 (as float)
      --lfsr-thres 0.01 (as float)
                        Retained for CLI compatibility; lfsr is now a per-
                        variant column (no split/rejoin).
  export_top_loci_3:
    Workflow Options:
      --qtl-type VAL (as str, required)
                        Combine per-region .top_loci.bed.gz files into one
                        study-level file. Public-facing name:
                        {name}.{qtl_type}[.{variant_tag}].top_loci.bed.gz
      --variant-tag ''
      --combine auto
                        'auto' (default): skip step 3 if region_name or
                        region_list is set; 'yes': force run; 'no': force skip.
      --region-list . (as path)
                        Re-declared so this step can detect step-1 subset
                        filters.
      --region-name  (as list)
  overlap_qtl_gwas_1:
    Workflow Options:
      --per-chunk 100 (as int)
      --gwas-meta-path . (as path)
      --qtl-meta-path . (as path)
      --gwas-file-path ''
      --qtl-file-path ''
      --region-list . (as path)
                        Optional: focus on a region list (last column = region
                        id) and/or region names.
      --region-name  (as list)
      --signal-cutoff 0.0 (as float)
                        PIP cutoff for the top-loci overlap. 0 = intersect all
                        top-loci variants (the legacy behaviour, which had its
                        CS filter commented out); the coarse block pairing that
                        feeds coloc wants breadth, not just the strongest
                        signal.
  overlap_qtl_gwas_2:
    Workflow Options:
      --[no-]step1-only (default to False)
                        set it as True if you only want run first step (running
                        parallel) and combine them manually after all finished

Workflow implementation

The global parameters and the workflow step definitions.

[global]
# A region list file documenting the chr_pos_ref_alt of a susie_object
parameter: cwd = path("output")
parameter: study = str
name = study

## Path to work directory where output locates
## Containers that contains the necessary packages
parameter: container = ""
import re
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
parameter: entrypoint= ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "96h"
# Memory expected
parameter: mem = "6G"
# Number of threads
parameter: numThreads = 2
parameter: windows = 1000000
# use this function to edit memory string for PLINK input
from sos.utils import expand_size

[susie_to_tsv_1]
# Migrated: consumes a FineMappingResult RDS (fine_mapping.R output) instead of a
# raw susieR object. The three legacy outputs map to fine_mapping_export.R views:
#   variant.tsv = topLoci   (per-variant: pip, posterior mean/sd, logBF, cs, region)
#   lbf.tsv     = lbf       (wide per-variant x per-effect log Bayes factor matrix)
#   effect.tsv  = cs_summary (per credible set: size, purity, V, logBF, lead variant)
# The FineMappingResult carries its own region/trait, so `region_list` is accepted
# for CLI compatibility but is no longer used.
parameter: region_list = path()
parameter: rds_path = paths
input: rds_path, group_by = 1
output: f"{cwd}/{_input:bn}.variant.tsv", f"{cwd}/{_input:bn}.lbf.tsv", f"{cwd}/{_input:bn}.effect.tsv"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand = '${ }', stdout = f"{_output[0]:nn}.stdout", stderr = f"{_output[0]:nn}.stderr", container = container, entrypoint = entrypoint
    export_script=${modular_script_dir}/pecotmr_integration/fine_mapping_export.R
    Rscript $export_script --input ${_input} --view topLoci --signal-cutoff 0 --output ${_output[0]}
    Rscript $export_script --input ${_input} --view lbf --output ${_output[1]}
    Rscript $export_script --input ${_input} --view cs_summary --output ${_output[2]}

[sQTL_susie_to_tsv_1]
# Migrated: identical to susie_to_tsv but preserves the sQTL / leafcutter2
# filename sanitization ("*" -> "N"; ":"/"+" scrubbed from the task tag), since
# intron ids contain shell-hostile characters. Consumes a FineMappingResult RDS
# and writes the topLoci / lbf / cs_summary views.
parameter: region_list = path()
parameter: rds_path = paths
input: rds_path, group_by = 1
input_name = f"{_input:bn}".replace('*', 'N')
output: f"{cwd}/{input_name}.variant.tsv", f"{cwd}/{input_name}.lbf.tsv", f"{cwd}/{input_name}.effect.tsv"
tags = f'{step_name}_{_output[0]:bn}'.replace(':', '_').replace('+', 'ps')
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = tags
bash: expand = '${ }', stdout = f"{_output[0]:nn}.stdout", stderr = f"{_output[0]:nn}.stderr", container = container, entrypoint = entrypoint
    export_script=${modular_script_dir}/pecotmr_integration/fine_mapping_export.R
    Rscript $export_script --input ${_input} --view topLoci --signal-cutoff 0 --output ${_output[0]}
    Rscript $export_script --input ${_input} --view lbf --output ${_output[1]}
    Rscript $export_script --input ${_input} --view cs_summary --output ${_output[2]}

[fsusie_to_tsv_1]
# Migrated: consumes an fSuSiE FineMappingResult RDS.
#   variant.tsv = topLoci (per-variant: pip, cs, cs_95_purity)
#   lbf.tsv     = lbf     (fSuSiE per-variant x per-effect lBF matrix)
# (Fixes the legacy bug where lbf.tsv was written the variant table, not the lbf
# matrix.) is_dummy = cs_95_purity < 0.5; the functional effect curves + peak
# positions are exported by fsusie_extract_effect / fsusie_affected_region
# (credible_band / affected_regions views), so they are not duplicated here.
parameter: region_list = path()
parameter: rds_path = paths
input: rds_path, group_by = 1
output: f"{cwd}/{_input:bn}.variant.tsv", f"{cwd}/{_input:bn}.lbf.tsv"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output[0]:bn}"
bash: expand = "${ }", stdout = f"{_output[0]:nn}.stdout", stderr = f"{_output[0]:nn}.stderr", container = container, entrypoint = entrypoint
    export_script=${modular_script_dir}/pecotmr_integration/fine_mapping_export.R
    Rscript $export_script --input ${_input} --view topLoci --signal-cutoff 0 --output ${_output[0]}
    Rscript $export_script --input ${_input} --view lbf --output ${_output[1]}

[*_to_tsv_2]
# Reduce: concatenate the per-region variant tables (the *variant.tsv outputs of
# the parallel _1 map step) into one file, in R via data.table::rbindlist
# (fill=TRUE aligns columns). _input carries all of _1's outputs, so keep only
# the variant tables (not the lbf / effect ones).
parameter: name = f'{_input[0]:b}'.split(".")[0]
input: group_by = "all"
output: f"{cwd}/{name}.all_variants.tsv"
R: expand = "${ }", stdout = f"{_output:n}.stdout", stderr = f"{_output:n}.stderr", container = container, entrypoint = entrypoint
    library(data.table)
    files <- Filter(function(f) endsWith(f, "variant.tsv"), c(${_input:r,}))
    fwrite(rbindlist(lapply(files, fread), fill = TRUE), ${_output:r}, sep = "\t")

[susie_tsv_collapse]
# Concatenate the per-region lbf tables (named *.chr<N>_<start>_<end>...lbf.tsv)
# by chromosome. Chromosome split is plain python (no pandas); the concat is R
# data.table::rbindlist (headers handled by fread; fill=TRUE aligns columns).
parameter: tsv_path = paths
chromosome = list(set(str(x).split(".")[-5].split("_")[0].replace("chr", "") for x in tsv_path))
input: tsv_path, for_each = "chromosome"
output: f'{cwd}/{_input[0]:bnnnnnnn}.chr{_chromosome}.unisusie_rss.lbf.tsv'
R: expand = "${ }", stdout = f"{_output}.stdout", stderr = f"{_output}.stderr", container = container, entrypoint = entrypoint
    library(data.table)
    files <- list.files(dirname("${_input[0]}"),
                        pattern = "\\.chr${_chromosome}_.*lbf\\.tsv$", full.names = TRUE)
    fwrite(rbindlist(lapply(files, fread), fill = TRUE), "${_output}", sep = "\t")

[susie_pip_landscape_plot]
# Per-region PIP landscape from a FineMappingResult (fine_mapping output).
# Modernized from the legacy multi-TSV plot: pecotmr does no plotting, so the
# ggplot lives in fine_mapping_pip_plot.R. (The regulatory/gene annotation
# overlay is a deferred enhancement; annot_tibble is kept for CLI compat.)
parameter: rds_path = paths
parameter: annot_tibble = path("~/Annotatr_builtin_annotation_tibble.tsv")
input: rds_path, group_by = 1
output: f'{cwd}/{_input:bn}.pip_landscape_plot.pdf'
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, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_pip_plot.R \
        --input ${_input} \
        --output ${_output}

[susie_upsetR_plot]
# UpSet plot of credible-set variant overlap across the FineMappingResult(s)'
# (study, context, trait, method) tuples. pecotmr does no plotting; the UpSet
# lives in fine_mapping_upset.R.
parameter: rds_path = paths
input: rds_path
output: f'{cwd}/{name}.UpSetR.pdf'
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, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_upset.R \
        --input ${" ".join([str(x) for x in _input])} \
        --output ${_output}

[susie_upsetR_cs_plot]
# Credible-set overlap UpSet plot across the FineMappingResult(s), via
# fine_mapping_upset.R (pecotmr does no plotting). trait_to_select is retained
# for CLI compatibility; the single-trait "CS shared with other phenotypes"
# sub-view is a deferred wrapper option.
parameter: rds_path = paths
parameter: trait_to_select = 1
input: rds_path
output: f'{cwd}/{name}.UpSetR_cs.pdf'
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, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_upset.R \
        --input ${" ".join([str(x) for x in _input])} \
        --output ${_output}

[tmp_annotatation_of_snps_1]
parameter: SNP_list = path
parameter: annotation_rds = path
input: SNP_list
output: f'{cwd}/{_input:b}.annotated.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
R: expand = '${ }', stdout = f"{_output}.stdout", stderr = f"{_output}.stderr", container = container
    library("dplyr")
    library("readr") 
    library("purrr")
    library("stringr")
    sharing_snp = readRDS("${_input}")
    sharing_snp_fsusie = sharing_snp[[1]]%>%filter(haQTL == 1 | mQTL == 1)
    sharing_snp_fsusie = sharing_snp_fsusie%>%mutate(X1 =  read.table(text = sharing_snp_fsusie$variants, sep = ":")$V1, X2 = read.table(text = read.table(text = sharing_snp_fsusie$variants, sep = ":")$V2 , sep = "_")$V1  )
    sharing_snp_fsusie = sharing_snp_fsusie%>%select(variants,chr = X1, pos = X2)
    annotation = readRDS("${annotation_rds}")
    print("data loaded")
    result = sharing_snp_fsusie%>%mutate(annot = map2( chr,pos , ~ annotation%>%filter(X1 == .x, X2 <= .y, X3 >= .y)%>%pull(X5)))%>%mutate(annot = map_chr(annot, ~paste0(.x ,collapse = ",")) )
    print("snp annotated")
    result%>%saveRDS("${_output}")

[tmp_annotatation_of_snps_2]
parameter: SNP_list = path
parameter: annotation_rds = path
output: f'{cwd}/{_input:b}.annotated_rev.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
R: expand = '${ }', stdout = f"{_output}.stdout", stderr = f"{_output}.stderr", container = container
    library("dplyr")
    library("readr") 
    library("purrr")
    library("stringr")
    result = readRDS(${_input:r})
    result_rev = tibble(annot = unique(annotation$X5))%>%mutate(variants = map(annot, ~  result%>%filter( str_detect(annot,.x))%>%pull(variants)) )%>%mutate( variants = map_chr(variants,~paste0(.x ,collapse = ","))  )
    result_rev%>%saveRDS("${_output}")

[fsusie_extract_effect]
# Fitted effect curve + credible band per fSuSiE credible set, via
# fsusieCredibleBand (fine_mapping_export.R `credible_band` view). pecotmr does
# the wavelet credible band; it needs an UNtrimmed fSuSiE fit (fine_mapping
# trim=FALSE). The effect-curve PDF is a deferred plot wrapper (pecotmr does no
# plotting); annot_tibble is kept for CLI compatibility.
parameter: rds_path = paths
parameter: annot_tibble = path("~/Annotatr_builtin_annotation_tibble.tsv")
input: rds_path, group_by = 1
output: f'{cwd}/{_input:bn}.estimated_effect.tsv'
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, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_export.R \
        --input ${_input} \
        --view credible_band \
        --output ${_output}

[fsusie_affected_region]
# GRanges of intervals where the fSuSiE credible band excludes 0 (with a
# direction mcol), via fsusieAffectedRegions (fine_mapping_export.R
# `affected_regions` view). Needs an UNtrimmed fSuSiE fit (fine_mapping
# trim=FALSE). The region PDF is a deferred plot wrapper.
parameter: rds_path = paths
input: rds_path, group_by = 1
output: f'{cwd}/{_input:bn}.affected_region.tsv'
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, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_export.R \
        --input ${_input} \
        --view affected_regions \
        --output ${_output}

#[uni_susie_2]
input: group_with = "genoFile"
output: f"{_input:n}.vcf.bgz"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
R: expand = '${ }', stdout = f"{_output:nn}.stdout", stderr = f"{_output:nn}.stderr", container = container
   ## Define create_vcf function
           create_vcf = function (chrom, pos, nea, ea, snp = NULL, ea_af = NULL, effect = NULL, 
        se = NULL, pval = NULL, name = NULL,cs = NULL, pip = NULL) 
    {
        stopifnot(length(chrom) == length(pos))
        if (is.null(snp)) {
            snp <- paste0(chrom, ":", pos)
        }
        snp <- paste0(chrom, ":", pos)
        nsnp <- length(chrom)
        gen <- list()
        ## Setupt data content for each sample column
        if (!is.null(ea_af)) 
            gen[["AF"]] <- matrix(ea_af, nsnp)
        if (!is.null(effect)) 
            gen[["ES"]] <- matrix(effect, nsnp)
        if (!is.null(se)) 
            gen[["SE"]] <- matrix(se, nsnp)
        if (!is.null(pval)) 
            gen[["LP"]] <- matrix(-log10(pval), nsnp)
        if (!is.null(cs)) 
            gen[["CS"]] <- matrix(cs, nsnp)
        if (!is.null(pip)) 
            gen[["PIP"]] <- matrix(pip, nsnp)
        gen <- S4Vectors::SimpleList(gen)
        
      ## Setup snps info for the fix columns
        gr <- GenomicRanges::GRanges(chrom, IRanges::IRanges(start = pos, 
            end = pos + pmax(nchar(nea), nchar(ea)) - 1, names = snp))
         coldata <- S4Vectors::DataFrame(Studies = name, row.names = name)
    ## Setup header informations
        hdr <- VariantAnnotation::VCFHeader(header = IRanges::DataFrameList(fileformat = S4Vectors::DataFrame(Value = "VCFv4.2", 
            row.names = "fileformat")), sample = name)
        VariantAnnotation::geno(hdr) <- S4Vectors::DataFrame(Number = c("A", 
            "A", "A", "A", "A", "A"), Type = c("Float", "Float", 
            "Float", "Float", "Float", "Float"), Description = c("Effect size estimate relative to the alternative allele", 
            "Standard error of effect size estimate", "-log10 p-value for effect estimate",  
            "Alternate allele frequency in the association study",
            "The CS this variate are captured, 0 indicates not in any cs", "The posterior inclusion probability to a CS"), 
            row.names = c("ES", "SE", "LP", "AF", "CS", "PIP"))
      ## Save only the meta information in the sample columns 
        VariantAnnotation::geno(hdr) <- subset(VariantAnnotation::geno(hdr), 
            rownames(VariantAnnotation::geno(hdr)) %in% names(gen))
      ## Save VCF 
        vcf <- VariantAnnotation::VCF(rowRanges = gr, colData = coldata, 
            exptData = list(header = hdr), geno = gen)
        VariantAnnotation::alt(vcf) <- Biostrings::DNAStringSetList(as.list(ea))
        VariantAnnotation::ref(vcf) <- Biostrings::DNAStringSet(nea)
      ## Add fixed values
        VariantAnnotation::fixed(vcf)$FILTER <- "PASS"
          return(sort(vcf))
        }
    library("susieR")
    library("dplyr")
    library("tibble")
    library("purrr")
    library("readr")
    library("tidyr")
    library("stringr")
    
    # Get list of cs snps
    susie_list = readRDS(${_input:r})
    susie_tb_ls = list()
    for (i in 1:length(susie_list)){
        susie_tb = tibble( snps =  names(susie_list[[1]]$pip)[which( susie_list[[i]]$pip >= 0)], snps_index = which(( susie_list[[i]]$pip >= 0))  )
        susie_tb_ls[[i]]= susie_tb%>%mutate( cs = map(snps_index,~which( susie_list[[i]]$sets$cs %in% .x))%>%as.numeric%>%replace_na(0),
                                 pip = map_dbl(snps_index,~( susie_list[[i]]$pip[.x])),
                                 coef = map_dbl(snps_index,~(coef.susie( susie_list[[i]])[.x+1])))
        }
    if(length(susie_tb_ls) >= 2){ 
      for(i in 2:length(susie_tb_ls)){
          susie_tb_ls[[i]] = full_join(susie_tb_ls[[i-1]],susie_tb_ls[[i]], by = "snps") 
        }
    }
    m = c("cs","pip","coef")    
    output = list()
    for(i in m){
    output[[i]] = susie_tb_ls[[length(susie_tb_ls)]]%>%select(contains(i))%>%as.matrix
    }
    snps_tb = susie_tb_ls[[length(susie_tb_ls)]]%>%mutate(
                         chr = map_chr(snps,~read.table(text = .x,sep = ":",as.is = T)$V1),
                         pos_alt_ref = map_chr(snps,~read.table(text = .x,sep = ":",as.is = TRUE)$V2),
                         pos = map_dbl(pos_alt_ref,~read.table(text = .x,sep = "_",as.is = TRUE)$V1),
                         alt = map_chr(pos_alt_ref,~read.table(text = .x,sep = "_",as.is = TRUE, colClass = "character")$V2),
                         ref = map_chr(pos_alt_ref,~read.table(text = .x,sep = "_",as.is = TRUE, colClass = "character")$V3))
    
    snps_tb = snps_tb%>%filter(str_detect(ref, "[ACTG]") & str_detect(alt, "[ACTG]"))
    output_vcf = create_vcf(
            chrom = snps_tb$chr,
             pos = snps_tb$pos,
             ea = snps_tb$alt,
             nea = snps_tb$ref,
             effect = snps_tb%>%select(contains("coef"))%>%as.matrix ,
             pip = snps_tb%>%select(contains("pip"))%>%as.matrix,
             cs = snps_tb%>%select(contains("cs"))%>%as.matrix,
             name = names(susie_list))
    VariantAnnotation::writeVcf(output_vcf,${_output:nr},index = TRUE)

[mv_susie_2]
# Write per-context fine-mapping VCFs (ES / SE / LP / AF + PIP / CS / within_cs_pip
# in the sample column) from a FineMappingResult via fine_mapping_vcf.R ->
# pecotmr::writeSumstatsVcf, replacing the inline create_vcf. writeSumstatsVcf
# splits a multi-context collection into one .vcf.bgz per context; the runtime
# set is tracked with output: dynamic() (so {cwd} must be created up front, since
# a dynamic output -- unlike a concrete target -- does not pre-create its dir).
parameter: rds_path = paths
import os
os.makedirs(str(cwd), exist_ok=True)
input: rds_path, group_by = 1
output: dynamic(f"{cwd}/{_input:bn}.*.vcf.bgz")
task: trunk_workers = 1, trunk_size = 1, walltime = '2h', mem = mem, cores = 1, tags = f'{step_name}_{_input:bn}'
bash: expand = '${ }', stdout = f"{cwd}/{_input:bn}.vcf.stdout", stderr = f"{cwd}/{_input:bn}.vcf.stderr", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_vcf.R \
        --input ${_input} \
        --output ${cwd}/${_input:bn}.vcf.bgz \
        --split-by-context

# Exporting cis susie_twas results
[cis_results_export_1, gwas_results_export_1, export_top_loci_1]
# per chunk we process at most 200 datasets
parameter: per_chunk = 200
# Region list should have last column being region name.
parameter:region_file=path()
# the path stored original output files
parameter:file_path=''
# assuming orinal files are named as prefix.id.suffix
# prefix of the original files (before id) to identify file.
parameter:prefix=[]
# suffix of the original files (after id) to identify file.
parameter:suffix=str
# if need to rename context, please load condition meta file
parameter: condition_meta = path()
# if the pip all variants in a cs < this threshold, then remove this cs
parameter: pip_thres = 0.05
# only keep top cs_size variants in one cs
parameter: cs_size = 3
# provide exported meta to filter the exported genes
parameter: exported_file = 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
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 = []
# the (aligned) geno bim file to check allele flipping. Retained for CLI
# compatibility; allele alignment is now done upstream by fine_mapping.R
# (matchVariants), so it is unused here.
parameter: geno_ref = path()
# context meta file to map renamed context backs
parameter: context_meta = path()
# CS purity (min.abs.corr) cutoff(s); the first value is used for conditions_top_loci.
parameter: min_corr = []
# set this parameter as `True` when exporting gwas data.
parameter: gwas = False
# set this parameter as `True` when exporting fsusie data.
parameter: fsusie = False
# set this parameter as `True` when exporting metaQTL data.
parameter: metaQTL = False
# set this parameter as `True` when exporting mnm data.
parameter: mnm = False

# Orchestration only (no pandas / no helper defs): resolve, per region, the
# existing per-study/context FineMappingResult input file(s). The heavy legacy
# work (allele alignment, top-loci processing, combined-db assembly) now lives
# in fine_mapping.R + the fine_mapping_cis_db_export.R wrapper below.
import os

_region_rows = [tuple((_l.rstrip('\n').split('\t') + ['', '', '', ''])[:4])
                for _l in open(region_file)
                if _l.strip() and not _l.startswith('#')]

_focus = set(region_name)
if region_list.is_file():
    _focus |= {_l.rstrip('\n').split('\t')[-1] for _l in open(region_list)
               if _l.strip() and not _l.startswith('#')}
if len(_focus) > 0:
    _region_rows = [r for r in _region_rows if r[3] in _focus or f"{r[0]}_{r[3]}" in _focus]

# candidate names: prefix.id.suffix (plain id first, then chr_id), pick existing.
_pfxs = list(prefix) if len(prefix) > 0 else ['']
_meta = []
_data = []
for _chr, _start, _end, _id in _region_rows:
    _plain = [f"{file_path}/{(p + '.') if p else ''}{_id}.{suffix}" for p in _pfxs]
    _plain = [c for c in _plain if os.path.exists(c)]
    _chrid = [f"{file_path}/{(p + '.') if p else ''}{_chr}_{_id}.{suffix}" for p in _pfxs]
    _chrid = [c for c in _chrid if os.path.exists(c)]
    _files = _plain if _plain else _chrid
    if _files:
        _meta.append((_chr, _start, _end, _id))
        _data.append(_files)

regional_data = {'meta': _meta, 'data': _data}
meta_info = regional_data['meta']
stop_if(len(regional_data['data']) == 0, 'No input FineMappingResult files found for any region.')
# db suffix mirrors cis_results_export_2's block_results_db / cis_results_db split
_db_kind = 'block_results_db' if gwas else 'cis_results_db'
_purity_arg = f"--min-purity {min_corr[0]}" if len(min_corr) > 0 else ""

input: regional_data["data"], group_by = lambda x: regional_data["data"], group_with = "meta_info"
output: f"{cwd}/{name}_cache/{name}_{_meta_info[3]}.tsv"
task: trunk_workers = job_size, walltime = walltime, trunk_size = job_size, mem = mem, cores = numThreads, tags = f'{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_cis_db_export.R \
        --input ${_input} \
        --combined-data-output ${cwd}/${name}.${_meta_info[3]}.${_db_kind}.rds \
        ${('--combined-data-sumstats-output %s/%s.%s.%s.sumstats.rds' % (cwd, name, _meta_info[3], _db_kind)) if not (fsusie or mnm) else ''} \
        --region-id ${_meta_info[3]} --chr ${_meta_info[0]} --start ${_meta_info[1]} --end ${_meta_info[2]} \
        --signal-cutoff ${pip_thres} ${_purity_arg} \
        --meta-output ${_output} \
        --pip-sum-output ${cwd}/${name}.${_meta_info[3]}.pip_sum

[cis_results_export_2, gwas_results_export_2]
# set it as True if you only want run first step (running parallel) and combine them manually after all finished 
parameter: step1_only = False
skip_if(step1_only)
# provide exported meta to filter the exported genes 
parameter: exported_file = path()
# optional: qtl or gwas, there is slightly different in qtl and gwas rds file
parameter: gwas = False
input: group_by = 'all'
output: f"{cwd}/{name}.{'block_results_db' if gwas else 'cis_results_db'}.tsv" 
# stop_if(_input[0] not in locals().keys(), 'All files have been exported already') #FIXME should we remove to a separate file. sothat we can stop globally as above
task: trunk_workers = 1, walltime = '1h', trunk_size = 1, mem = '16G', cores = 1, tags = f'{_output:bn}'
bash: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    if [ -e "${_output:ad}/${name}_cache/" ]; then
        sed 's/^chr/#chr/' `ls ${_output:ad}/${name}_cache/*tsv |head -n1` | head -n 1 > ${_output:an}.temp
        tail -n +2 -q ${_output:ad}/${name}_cache/*.tsv >> ${_output:an}.temp
        error_files=$(find "${_output:ad}/${name}_cache/" -type f -name "*_error")

        if [[ -n $error_files ]]; then
            cat $error_files >> ${_output:an}.error_genes
        else
            echo "No truncated files detected"
        fi
    else
        echo "All files have been exported already"
    fi
    
R: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    if (file.exists(paste0(${_output:anr},".temp"))) {
        library(tidyverse)
        meta <- read_delim(paste0(${_output:anr},".temp"), delim = '\t')

        if (${"TRUE" if exported_file.is_file() else "FALSE"}){
          exp_meta <- read_delim(${exported_file:r}, delim = '\t')
          meta <- bind_rows(meta, exp_meta) %>%
              group_by(`#chr`, start, end, region_id, TSS) %>%
              summarise(across(c(original_data, combined_data, combined_data_sumstats, conditions, conditions_top_loci), 
                               ~paste(unique(.), collapse = ",")),
                        .groups = 'drop')
              }

        write_delim(meta, ${_output:r}, delim = '\t')
    }

[cis_results_export_3]
# set it as True if you only want run first step (running parallel) and combine them manually after all finished 
parameter: step1_only = False

skip_if(step1_only)
bash: expand = "${ }", container = container, stderr = f'{_input:n}.stderr', stdout = f'{_input:n}.stdout'   
    rm -rf "${_input:ad}/${name}_cache/"
    rm -rf ${_input:an}.temp

#get union of step1
#1200 blocks costed ~5mins with in one for loop
[gwas_results_export_3]
parameter: step1_only = False

skip_if(step1_only)
output: f"{cwd}/{name}.union_export.tsv.gz" 
task: trunk_workers = 1, walltime = '1h', trunk_size = 1, mem = '16G', cores = 1, tags = f'{_output:bn}'
bash: expand = "${ }", container = container, stderr = f'{_input:n}.stderr', stdout = f'{_input:n}.stdout'   
    rm -rf "${_input:ad}/${name}_cache/"
    rm -rf ${_input:an}.temp

R: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    library(tidyverse)
    library(data.table)

    mtx <- read_delim(${_input:r})
    files <- mtx %>% filter(!is.na(conditions_top_loci)) %>% pull(combined_data) %>% paste0(${_input[0]:dr},'/',.)
    all_top_loci <- data.frame()

    for (i in seq_along(files)) {
      file <- files[i]
      res <- readRDS(file)

      file_top_loci <- lapply(names(res), function(block) {
        lapply(names(res[[block]]), function(study) {
          lapply(names(res[[block]][[study]]), function(method) {
            if (!is.null(res[[block]][[study]][[method]]$top_loci)) {
              temp_df <- res[[block]][[study]][[method]]$top_loci
              mutate(temp_df, study = study, method = method, block = block)
            } else {
              NULL  
            }
          })
        }) %>% bind_rows()  
      }) %>% bind_rows() 

      all_top_loci <- bind_rows(all_top_loci, file_top_loci)
      if (i %% 100 == 0) {
        message(sprintf("Have processed %d files.", i))
      }
    }
    fwrite(all_top_loci, ${_output:r})

# simply combine seperate meta files, not works for having an exsiting exported file for now. 
[combine_export_meta]
parameter: cache_path=path
parameter: output_file = str
parameter: remove_cache = False
output: f"{cache_path:d}/{output_file}" 
bash: expand = "${ }", container = container, stderr = f'{_input:n}.stderr', stdout = f'{_input:n}.stdout'   
    head -n 1 -q ${cache_path}/*.tsv | sed 's/^chr/#chr/' | head -n 1 > ${_output}
    tail -n +2 -q ${cache_path}/*.tsv >> ${_output}
    if ${"true" if remove_cache else "false"}; then
        rm -rf ${cache_path}
    fi

[export_top_loci_2]
# Export the per-region top-loci BED from each cis_results_db (an FMR in the
# modern pipeline) via the pecotmr-backed fine_mapping_top_loci_bed.R wrapper.
# getTopLoci does the PIP/CS + independent purity filtering; lfsr /
# conditional_effect are per-(variant, context) numeric columns, so the legacy
# filter_lfsr semicolon split/rejoin is gone.
# `export_suffix` deliberately differs from step 1's `--suffix` (input rds) to avoid CLI collision.
parameter: export_path = cwd
parameter: export_prefix = ''  # falls back to step-1 `name` if empty
parameter: export_suffix = ''  # empty => derived from `gwas` below
# Re-declared so the suffix can default by mode: step 1 writes
# block_results_db.rds under --gwas True, else cis_results_db.rds.
parameter: gwas = False
parameter: fsusie_prefix = ''
parameter: preset_top_loci = False
# PIP cutoff for getTopLoci; optional independent CS purity (min.abs.corr) cutoff (<0 => none).
parameter: signal_cutoff = 0.025
parameter: min_purity = -1.0
# Retained for CLI compatibility; lfsr is now a per-variant column (no split/rejoin).
parameter: lfsr_thres = 0.01

_pfx = export_prefix if export_prefix else name
import glob as _glob
_suffix = export_suffix if export_suffix else ('block_results_db.rds' if gwas else 'cis_results_db.rds')
_rds_files = sorted(_glob.glob(f"{export_path}/{_pfx}.{fsusie_prefix}*.{_suffix}"))
stop_if(len(_rds_files) == 0, f"No RDS files matched: {export_path}/{_pfx}.{fsusie_prefix}*.{_suffix}")
_purity_arg = f"--min-purity {min_purity}" if float(min_purity) >= 0 else ""

input: _rds_files, group_by = 1
output: f"{cwd}/summary/{_input:bnn}.top_loci.bed.gz"
task: trunk_workers = 1, walltime = '1h', trunk_size = job_size, mem = '16G', cores = 1, tags = f'{_output:bn}'
bash: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', entrypoint=entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_top_loci_bed.R \
        --input ${_input} \
        --signal-cutoff ${signal_cutoff} ${_purity_arg} \
        --output ${_output}

[export_top_loci_3]
# Combine per-region .top_loci.bed.gz files into one study-level file.
# Public-facing name: {name}.{qtl_type}[.{variant_tag}].top_loci.bed.gz
parameter: qtl_type = str
parameter: variant_tag = ''
# 'auto' (default): skip step 3 if region_name or region_list is set; 'yes': force run; 'no': force skip.
parameter: combine = 'auto'
# Re-declared so this step can detect step-1 subset filters.
parameter: region_list = path()
parameter: region_name = []
stop_if(combine not in ('auto', 'yes', 'no'), f"--combine must be auto|yes|no, got {combine}")
_region_filter_set = (region_list.is_file() or len(region_name) > 0)
_do_combine = (combine == 'yes') or (combine == 'auto' and not _region_filter_set)
skip_if(not _do_combine, f"Skipping combine: combine={combine}, region_filter_set={_region_filter_set}")
_tag = f".{variant_tag}" if variant_tag else ""
input: group_by = 'all'
output: f"{cwd}/summary/{name}.{qtl_type}{_tag}.top_loci.bed.gz"
task: trunk_workers = 1, walltime = '2h', trunk_size = 1, mem = '8G', cores = 1, tags = f'{_output:bn}'
bash: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', entrypoint=entrypoint
    set -euo pipefail
    mkdir -p "$(dirname "${_output}")"
    GZ=$(command -v bgzip || echo gzip)
    {
        # header from first non-empty file
        for f in ${_input}; do
            if [ "$(zcat -f "$f" 2>/dev/null | wc -l)" -ge 2 ]; then
                zcat -f "$f" | head -1
                break
            fi
        done
        # data rows from every non-empty file
        for f in ${_input}; do
            if [ "$(zcat -f "$f" 2>/dev/null | wc -l)" -ge 2 ]; then
                zcat -f "$f" | tail -n +2
            fi
        done
    } | $GZ > ${_output}

[overlap_qtl_gwas_1]
parameter: per_chunk = 100
parameter: gwas_meta_path = path()
parameter: qtl_meta_path = path()
parameter: gwas_file_path = ''
parameter: qtl_file_path = ''
# Optional: focus on a region list (last column = region id) and/or region names.
parameter: region_list = path()
parameter: region_name = []
# PIP cutoff for the top-loci overlap. 0 = intersect all top-loci variants (the
# legacy behaviour, which had its CS filter commented out); the coarse block
# pairing that feeds coloc wants breadth, not just the strongest signal.
parameter: signal_cutoff = 0.0
# Orchestration only (no pandas / no helper defs): pair each QTL region with the
# GWAS blocks it overlaps by coordinate, resolving both sides' original_data
# (per-study/per-block S4 FineMappingResult files). fine_mapping_overlap.R does
# the allele-aware variant overlap + stamps block_top_loci onto the meta row.
import os, csv

if qtl_file_path == '':  qtl_file_path  = str(qtl_meta_path.parent)
if gwas_file_path == '': gwas_file_path = str(gwas_meta_path.parent)

_hasTop = lambda r: (r.get('conditions_top_loci') or '').strip() not in ('', 'NA')
_g = [r for r in csv.DictReader(open(gwas_meta_path), delimiter='\t') if _hasTop(r)]
_q = [r for r in csv.DictReader(open(qtl_meta_path),  delimiter='\t') if _hasTop(r)]

_focus = set(region_name)
if region_list.is_file():
    _focus |= {l.rstrip('\n').split('\t')[-1] for l in open(region_list) if l.strip() and not l.startswith('#')}
if len(_focus) > 0:
    _q = [r for r in _q if r['region_id'] in _focus or f"{r['#chr']}_{r['region_id']}" in _focus]

_meta = []
_qtl_data = []
for q in _q:
    qchr, qs, qe = q['#chr'], int(q['start']), int(q['end'])
    blocks = [g for g in _g if g['#chr'] == qchr and int(g['start']) <= qe and int(g['end']) >= qs]
    gfiles = [f"{gwas_file_path}/{x.strip()}" for g in blocks for x in g['original_data'].split(',') if x.strip()]
    gfiles = [f for f in gfiles if os.path.exists(f)]
    qfiles = [f"{qtl_file_path}/{x.strip()}" for x in q['original_data'].split(',') if x.strip()]
    qfiles = [f for f in qfiles if os.path.exists(f)]
    if gfiles and qfiles:
        _meta.append((qchr, q['start'], q['end'], q['region_id'], gfiles))
        _qtl_data.append(qfiles)

stop_if(len(_qtl_data) == 0, 'No overlapping QTL x GWAS region pairs with existing files.')
meta_info = _meta

input: [f for grp in _qtl_data for f in grp], group_by = lambda x: _qtl_data, group_with = "meta_info"
output: f"{cwd}/gwas_qtl/cache/{name}_gwas_batch_meta_{_meta_info[3]}.tsv"
task: trunk_workers = 1, walltime = '1h', trunk_size = job_size, mem = '16G', cores = 1, tags = f'{_output:bn}'
bash: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_overlap.R \
        --qtl ${" ".join([str(x) for x in _input])} \
        --gwas ${" ".join(_meta_info[4])} \
        --signal-cutoff ${signal_cutoff} \
        --qtl-meta ${qtl_meta_path} --region-id ${_meta_info[3]} \
        --meta-output ${_output} \
        --output ${cwd}/gwas_qtl/cache/${name}_gwas_batch_export_${_meta_info[3]}.tsv

[overlap_qtl_gwas_2]
# set it as True if you only want run first step (running parallel) and combine them manually after all finished 
parameter: step1_only = False
skip_if(step1_only)
input: group_by = 'all'
output: f"{cwd}/{name}.overlapped.gwas.tsv"
task: trunk_workers = 1, walltime = '1h', trunk_size = job_size, mem = '16G', cores = 1, tags = f'{_output:bn}'
R: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    library(data.table)
    exp_path <- ${_input[0]:adr}
    meta_files <- c(${",".join(['"%s"' % x.absolute() for x in _input])})
    exp_files <- list.files(exp_path, "_gwas_batch_export", full.names = T)
    meta_list <- exp_list <- list()
    meta_combined <- rbindlist(lapply(meta_files, fread), fill = TRUE)
    exp_combined <- rbindlist(lapply(exp_files, fread), fill = TRUE)
    fwrite(exp_combined, gsub("tsv","export.csv.gz","${_output}"))
    fwrite(meta_combined, "${_output}", sep = '\t')
  
#bash: expand = "${ }", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
    # rm -rf ${_input[0]:adr}