Fine-mapping result post-processing

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

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

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

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

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

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

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

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

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

sos run pipeline/mnm_postprocessing.ipynb -h

[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}