MASH analysis pipeline with posterior computation

MASH analysis pipeline with posterior computation#

Applies a fitted MASH model to compute posterior quantities for gene-SNP effect-size chunks and derives contrast and feature-score summaries.

Overview#

For each input chunk (a list of matrices bhat/sbhat/Z), the posterior workflow loads the MASH model and calls mash_compute_posterior_matrices. Additional workflows compute posterior contrasts between conditions and feature-level scores (meta, fine-mapped, n-significant, and p-value pairs) from the contrast results.

Fitting the mixture model and applying it are separate jobs. mash_fit learns which patterns of sharing exist across conditions; this notebook applies that fitted model to each chunk of effects, shrinking noisy estimates towards the patterns the data support. Effects that look condition-specific because of noise get pulled towards the shared pattern, and genuinely specific ones do not.

It is designed to run posterior computation in parallel over many input chunks.

When to run it. After mash_fit, on the same effect estimates. Posterior quantities - not the raw per-condition estimates - are what downstream comparisons should use.

Requires the MASH model from mash_fit (protocol_example.mash_model.rds).

Input#

posterior - computes per-region posterior mean/covariance:

  • --mash-model input/mash/protocol_example.mash_model.rds (a fitted MASH model RDS, the bare mash object, e.g. produced by mash_fit)

List of 9
 $ result           :List of 5
  ..$ PosteriorMean: num [1:2000, 1:34] 0.07633 -0.06986 0.06814 0.04333 0.15277 0.00605 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ PosteriorSD  : num [1:2000, 1:34] 0.28 0.437 0.304 0.299 0.42 0.329 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ NegativeProb : num [1:2000, 1:34] 0.415 0.555 0.42 0.412 0.359 0.467 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ lfsr         : num [1:2000, 1:34] 0.44 0.445 0.447 0.503 0.374 0.52 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ PosteriorCov : num [1:34, 1:34, 1:2000] 0.078574 0.00178 -0.001302 0.000344 -0.000536 0.010921 ...
  .. ..- attr(*, "dimnames")=List of 3
 $ loglik           : num -98517
 $ vloglik          : num [1:2000, 1] -48.2 -51.6 -49.4 -44.5 -53.6 -46.2 ...
 $ null_loglik      : num [1:2000] -49 -53.1 -50 -44 -54.8 -46.1 ...
 $ alt_loglik       : num [1:2000, 1] -48.2 -51.6 -49.3 -44.6 -53.5 -46.2 ...
 $ fitted_g         :List of 4
  ..$ pi          : Named num [1:1169] 0.0508 0 0 0 0 0 ...
  .. ..- attr(*, "names")= chr [1:1169] "null" "XtX.1" "tFLASH_default.1" "FLASH_default_1.1" "FLASH_default_2.1" ...
  ..$ Ulist       :List of 73
  .. ..$ XtX                         : num [1:34, 1:34] 0.1763 -0.0223 0.0796 -0.0213 0.0865 -0.0216 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ tFLASH_default              : num [1:34, 1:34] 1 -0.0025 0.33493 0.00075 0.3636 0.00286 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_1             : num [1:34, 1:34] 9.99e-05 -1.23e-07 -1.36e-10 -2.99e-07 6.28e-11 -1.68e-07 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_2             : num [1:34, 1:34] 1.49e-04 2.17e-08 4.35e-08 -3.44e-08 1.26e-07 9.20e-09 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_3             : num [1:34, 1:34] 6.12e-05 1.28e-09 3.24e-09 2.84e-09 6.26e-09 2.59e-09 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_4             : num [1:34, 1:34] 4.96e-05 -7.28e-10 -2.49e-10 -1.50e-08 7.31e-11 -1.70e-10 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_5             : num [1:34, 1:34] 1.00 -1.13e-07 2.75e-01 8.71e-03 4.58e-01 7.49e-04 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_6             : num [1:34, 1:34] 9.52e-05 1.82e-11 -1.35e-10 4.70e-10 -2.54e-11 9.50e-12 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_7             : num [1:34, 1:34] 4.68e-05 -4.41e-10 2.77e-10 -5.70e-10 -2.21e-11 2.13e-09 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_8             : num [1:34, 1:34] 0.022 0.1477 0.0211 0.0674 0.0203 0.0766 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_9             : num [1:34, 1:34] 1.86e-04 5.22e-14 -4.29e-08 4.64e-09 -2.18e-09 5.02e-14 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_10            : num [1:34, 1:34] 2.76e-05 -8.13e-10 1.48e-10 2.10e-10 2.85e-10 -5.92e-09 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_11            : num [1:34, 1:34] 1.21e-03 2.27e-08 -1.58e-03 -1.22e-04 -3.18e-02 3.68e-04 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_12            : num [1:34, 1:34] 1.49e-03 3.86e-06 1.48e-03 1.40e-06 1.26e-03 1.44e-05 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_13            : num [1:34, 1:34] 6.92e-04 8.42e-06 2.20e-04 -1.15e-05 1.24e-04 1.46e-05 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_14            : num [1:34, 1:34] 4.36e-03 6.44e-05 6.04e-03 2.85e-02 2.34e-03 6.49e-02 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_15            : num [1:34, 1:34] 2.73e-01 -2.14e-10 -2.92e-11 8.84e-10 2.79e-09 -6.92e-10 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_16            : num [1:34, 1:34] 0.000401 0.000185 0.000203 0.000398 0.000173 0.000186 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_17            : num [1:34, 1:34] 1.01e-03 -1.15e-07 1.09e-06 -1.24e-06 -5.17e-07 2.19e-07 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ FLASH_default_18            : num [1:34, 1:34] 6.41e-05 -1.99e-10 4.53e-08 3.19e-06 1.39e-09 4.34e-10 ...
  .. .. ..- attr(*, "dimnames")=List of 2
... (truncated)

  • --analysis-units input/finemapping/protocol_example.analysis_units.txt (a text file whose first column lists paths to posterior-input RDS chunks, one region per line)


The per-chunk files this manifest points at are prepared beforehand, one set per batch:

- `*.batch_*.yaml` - the gene-SNP pairs of interest for that batch, identified elsewhere (for example by fine-mapping analysis).
- `*.batch_*.rds` - the univariate summary statistics for those gene-SNPs, extracted from the YAML; this is what `posterior` consumes.
- `*.batch_*.stdout` - records SNPs present in the fine-mapping results but absent from the original `fastqtl` output.
<path>/protocol_example.posterior_input.rds region1

  • --posterior-vhat-files input/twas/protocol_example.posterior_vhat.rds (matching residual-variance (vhat) matrix RDS file(s))

 num [1:34, 1:34] 1 0 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
  ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...

Produces a manifest file (mash_output_list_all, tab-separated id path pairs) pointing to the per-region posterior RDS outputs.

mash_posterior_contrast - computes contrasts from the posterior output:

  • --posterior-file: a tab-separated manifest (id path) where each path is a per-region posterior RDS produced by the posterior step above.

  • --sum-file: a tab-separated manifest (id path) where each path is the corresponding raw summary-statistics RDS for that region, containing top-level bhat/sbhat matrices (the same data used to compute the posterior). The id values must match between --posterior-file and --sum-file.

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

  • --data-table-name bhat (name of the effect-size table to read from each input chunk)

  • --exclude-condition (conditions to drop before computing posteriors; empty by default)

Output#

posterior_1

  • cache/{name}.{chunk}.posterior.rds - one serialized posterior object per input chunk, loadable with readRDS(). The input is split into chunks so the work can be spread across cluster nodes. Each object holds PosteriorMean, PosteriorSD, NegativeProb, lfsr and PosteriorCov (effects x conditions).

List of 5
$ PosteriorMean: matrix [300, 34] 0.137296644006322 -0.311361510169076 0.317082792980895 -0.166211879988538 0.18908057204664 ...
$ PosteriorSD: matrix [300, 34] 0.344426275519538 0.449586894071048 0.431655774119752 0.469980878283076 0.417647229394468 ...
$ lfdr: matrix [300, 34] 0.0781635470570518 0.00282142909610905 0.0300484413895379 0.00250476599353006 0.0514063370711267 ...
$ NegativeProb: matrix [300, 34] 0.337086778060348 0.743597948338302 0.226498052642181 0.635857857046459 0.31287430116626 ...
$ lfsr: matrix [300, 34] 0.4152503251174 0.256402051661698 0.256546494031719 0.364142142953541 0.364280638237386 ...

posterior_2

  • {name}.{output_suffix}.posterior_list - a manifest of the per-chunk posterior files; the suffix comes from --output-suffix.

mash_posterior_contrast_1

  • contrast/{name}.{chunk}.posterior_contrast.rds - per-chunk posterior contrasts between conditions.

mash_posterior_contrast_2

  • {name}.posterior_sum.csv - the contrast summary table.

posterior_cntrast_plot

  • {name}.posterior_sum.png - the contrast plot.

feature_score_meta, feature_score_finemap, feature_score_nsig and feature_pval_pair

  • <step>/cache/{name}.featurescore{N}.rds - per-chunck score

  • {name}.{step_name}.feature_score_sum.csv - a summary table

Minimal Working Example#

Step 1. Compute MASH posteriors for each input chunk listed in the analysis-units file#

Timing: ~30 sec (on toy dataset)

sos run pipeline/mash_posterior.ipynb posterior \
    --cwd output/mash_posterior \
    --analysis-units input/finemapping/protocol_example.analysis_units.txt \
    --mash-model input/mash/protocol_example.mash_model.rds \
    --posterior-vhat-files input/twas/protocol_example.posterior_vhat.rds \
    --data-table-name strong \
    --exclude-condition 1 3

Step 2. Compute posterior contrasts between conditions for the sliced data#

Rscript -e 'p <- strsplit(readLines("input/finemapping/protocol_example.analysis_units.txt")[1], "\\s+")[[1]][1]; saveRDS(readRDS(p)$strong, "output/mash_posterior/protocol_example.region1_sumstats.rds")'
printf 'region1\t%s\n' "$PWD/output/mash_posterior/cache/protocol_example.posterior_input.posterior.rds" > output/mash_posterior/posterior_manifest.txt
printf 'region1\t%s\n' "$PWD/output/mash_posterior/protocol_example.region1_sumstats.rds" > output/mash_posterior/sum_manifest.txt

sos run pipeline/mash_posterior.ipynb mash_posterior_contrast \
    --cwd output/mash_posterior \
    --posterior-file output/mash_posterior/posterior_manifest.txt \
    --sum-file output/mash_posterior/sum_manifest.txt

Step 3. Plot the posterior contrast results#

Timing: TBD (on toy dataset)

sos run pipeline/mash_posterior.ipynb posterior_contrast_plot \
    --cwd output/mash_posterior \
    --analysis-units input/finemapping/protocol_example.analysis_units.txt

Step 4. Compute meta feature scores from the contrast results#

Timing: TBD (on toy dataset)

sos run pipeline/mash_posterior.ipynb feature_score_meta \
    --cwd output/mash_posterior \
    --analysis-units input/finemapping/protocol_example.analysis_units.txt \
    --posterior-file input/protocol_example.posterior.rds \
    --sum-file input/protocol_example.sumstats.rds

Step 5. Compute feature scores from contrast results using fine-mapped eQTL/pQTL#

Timing: TBD (on toy dataset)

sos run pipeline/mash_posterior.ipynb feature_score_finemap \
    --cwd output/mash_posterior \
    --analysis-units input/finemapping/protocol_example.analysis_units.txt \
    --posterior-file input/protocol_example.posterior.rds \
    --sum-file input/protocol_example.sumstats.rds

Step 6. Compute n-significant feature scores from contrast results#

Timing: TBD (on toy dataset)

sos run pipeline/mash_posterior.ipynb feature_score_nsig \
    --cwd output/mash_posterior \
    --analysis-units input/finemapping/protocol_example.analysis_units.txt \
    --posterior-file input/protocol_example.posterior.rds \
    --sum-file input/protocol_example.sumstats.rds

Step 7. Compute p-value-pair feature scores from contrast results#

Timing: TBD (on toy dataset)

sos run pipeline/mash_posterior.ipynb feature_pval_pair \
    --cwd output/mash_posterior \
    --analysis-units input/finemapping/protocol_example.analysis_units.txt \
    --posterior-file input/protocol_example.posterior.rds \
    --sum-file input/protocol_example.sumstats.rds

Command Interface#

sos run pipeline/mash_posterior.ipynb -h

usage: sos run pipeline/mash_posterior.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:
  posterior
  mash_posterior_contrast
  posterior_cntrast_plot
  feature_score_meta
  feature_score_finemap
  feature_score_nsig
  feature_pval_pair

Global Workflow Options:
  --cwd output (as path)
  --modular-script-dir code/script (as path)
  --name test
  --cells  (as list)
                        Conditions (contexts); order matters
  --group1  (as list)
                        Condition groups: replicate populations of one cell type
                        share contrast weight
  --group2  (as list)
  --group3  (as list)
  --job-size 1 (as int)
  --container ''
  --per-chunk 1 (as int)
                        Number of analysis units per job
  --output-prefix ''
  --output-suffix all
  --effect-model EE
                        Exchangable effect (EE) or exchangable z-scores (EZ)
  --vhat simple
                        Vhat estimate identifier (e.g. simple / identity / mle)
  --data fastqtl_to_mash_output/FastQTLSumStats.mash.rds (as path)
  --p-cut 1e-05 (as float)
                        Significance cutoff for the contrast summary / n-sig
                        feature score
  --walltime 1h
  --mem 16G
  --numThreads 1 (as int)

Sections
  posterior_1:          Compute MASH posteriors per analysis unit (one posterior
                        RDS per region)
    Workflow Options:
      --analysis-units VAL (as path, required)
                        File listing per-region data RDS paths (col 1); each
                        carries bhat/sbhat matrices
      --bhat-table-name bhat
                        Effect-size / standard-error list elements inside each
                        region RDS
      --shat-table-name sbhat
      --exclude-condition  (as list)
                        Conditions (columns) to exclude; names or 1-based
                        indices
  posterior_2:          Collect the per-region posterior paths into a single
                        list
  mash_posterior_contrast_1: Per-region posterior contrasts (deviation +
                        pairwise) via mashPosteriorContrast
    Workflow Options:
      --analysis-units VAL (as path, required)
                        File listing per-region data RDS paths (same units as
                        posterior_1)
      --orig-key bhat
                        Effect-size list element inside each region RDS (aligned
                        to the posterior)
      --grouping-recipe ''
                        Optional file of comma-separated condition groups (one
                        per line)
  mash_posterior_contrast_2: Summarize contrast significance across regions ->
                        CSV
  mash_posterior_contrast_3, posterior_cntrast_plot: Plot the contrast
                        significance summary as a symmetric heatmap
  feature_score_meta_1: Feature score (meta) from per-region contrast results
    Workflow Options:
      --contrast-units VAL (as path, required)
                        File listing per-region posterior-contrast RDS paths
      --meta-method REML
                        Random-effects estimator (metafor via pecotmr)
  feature_score_finemap_1: Feature score (finemap) using credible sets from
                        fine-mapping
    Workflow Options:
      --contrast-units VAL (as path, required)
      --fine-mapping VAL (as path, required)
                        Fine-mapping table RDS with cs_order / pip / variants
                        columns
      --conditions ''
                        Conditions to score (default: all contexts present in
                        the contrast)
  feature_score_finemap_2: Merge per-chunk finemap feature scores into one table
  feature_score_nsig_1: Feature score (nsig) from per-region contrast results
    Workflow Options:
      --contrast-units VAL (as path, required)
                        File listing per-region posterior-contrast RDS paths
      --p-cutoff 1e-05 (as float)
                        Significance cutoff for the n-significant ratio
  feature_pval_pair_1:  Feature score (pval_pair) from per-region contrast
                        results
    Workflow Options:
      --contrast-units VAL (as path, required)
                        File listing per-region posterior-contrast RDS paths
      --meta-method REML
      --se-cutoff 0.001 (as float)
                        SE floor for the pairwise meta-analysis
  feature_pval_pair_2, feature_score_meta_2, feature_score_nsig_2: Merge per-
                        chunk feature scores into one table (meta / nsig /
                        pval_pair)

Workflow implementation#

[global]
parameter: cwd = path('./output')
parameter: modular_script_dir = path('code/script')
parameter: name = 'test'
# Conditions (contexts); order matters
parameter: cells = []
# Condition groups: replicate populations of one cell type share contrast weight
parameter: group1 = []
parameter: group2 = []
parameter: group3 = []
parameter: job_size = 1
parameter: container = ''
# Number of analysis units per job
parameter: per_chunk = 1
parameter: output_prefix = ''
parameter: output_suffix = 'all'
# Exchangable effect (EE) or exchangable z-scores (EZ)
parameter: effect_model = 'EE'
# Vhat estimate identifier (e.g. simple / identity / mle)
parameter: vhat = 'simple'
parameter: data = path("fastqtl_to_mash_output/FastQTLSumStats.mash.rds")
# Significance cutoff for the contrast summary / n-sig feature score
parameter: p_cut = 0.00001
parameter: walltime = '1h'
parameter: mem = '16G'
parameter: numThreads = 1
data = data.absolute()
cwd = cwd.absolute()
if len(output_prefix) == 0:
    output_prefix = f"{data:bn}"
vhat_data = file_target(f"{cwd:a}/{output_prefix}.{effect_model}.V_{vhat}.rds")
mash_model = file_target(f"{cwd:a}/{output_prefix}.{effect_model}.V_{vhat}.mash_model.rds")

posterior#

take all the 13K genes, with slice_method = True, conditions missing from a unit have their rows and columns dropped from the prior model

# Compute MASH posteriors per analysis unit (one posterior RDS per region)
[posterior_1]
# File listing per-region data RDS paths (col 1); each carries bhat/sbhat matrices
parameter: analysis_units = path
# Effect-size / standard-error list elements inside each region RDS
parameter: bhat_table_name = 'bhat'
parameter: shat_table_name = 'sbhat'
# Conditions (columns) to exclude; names or 1-based indices
parameter: exclude_condition = []
posterior_input = [line.split()[0] for line in open(analysis_units).readlines() if line.strip() and not line.strip().startswith('#')]
input: posterior_input, group_by = 1
output: f"{cwd}/cache/{name}.{_input:bn}.posterior.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_posterior.R \
        --data ${_input} \
        --bhat-key ${bhat_table_name} \
        --shat-key ${shat_table_name} \
        --vhat-data ${vhat_data} \
        --mash-model ${mash_model} \
        --effect-model ${effect_model} \
        --exclude-condition "${','.join([str(x) for x in exclude_condition])}" \
        --output ${_output}

# Collect the per-region posterior paths into a single list
[posterior_2]
input: group_by = "all"
output: f"{cwd}/{name}.{output_suffix}.posterior_list"
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    for f in ${_input}; do echo "$f"; done > ${_output}

mash_posterior_contrast#

  • Add group in this module

    e.g. add MiGA-Microglia-scRNA data to this analysis to supplement the Mic cell count shortage, but the MiGA data comes from four sites with four datasets, and here I handled it in such a way that the different Mic’s were always in the same state during the deviation contrast analysis. If Mic is the one be compared, each Mic sample would be set as (n_populations - 1)/n_Mic instead of n_populations

# Per-region posterior contrasts (deviation + pairwise) via mashPosteriorContrast
[mash_posterior_contrast_1]
# File listing per-region data RDS paths (same units as posterior_1)
parameter: analysis_units = path
# Effect-size list element inside each region RDS (aligned to the posterior)
parameter: orig_key = 'bhat'
# Optional file of comma-separated condition groups (one per line)
parameter: grouping_recipe = ''
contrast_units = [line.split()[0] for line in open(analysis_units).readlines() if line.strip() and not line.strip().startswith('#')]
input: contrast_units, group_by = 1
output: f"{cwd}/contrast/{name}.{_input:bn}.posterior_contrast.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_posterior_contrast.R \
        --posterior ${cwd}/cache/${name}.${_input:bn}.posterior.rds \
        --orig-data ${_input} \
        --orig-key ${orig_key} \
        --cells ${','.join(cells)} \
        --group1 "${','.join(group1)}" \
        --group2 "${','.join(group2)}" \
        --group3 "${','.join(group3)}" \
        --grouping-recipe "${grouping_recipe}" \
        --output ${_output}

# Summarize contrast significance across regions -> CSV
[mash_posterior_contrast_2]
input: group_by = "all"
output: f"{cwd}/{name}.posterior_sum.csv"
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_posterior_contrast_summary.R \
        --contrast ${_input} \
        --cells ${','.join(cells)} \
        --p-cutoff ${p_cut} \
        --output ${_output}

# Plot the contrast significance summary as a symmetric heatmap
[mash_posterior_contrast_3, posterior_cntrast_plot]
input: group_by = "all"
output: f"{cwd}/{name}.posterior_sum.png"
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_posterior_contrast_plot.R \
        --data ${_input} \
        --output ${_output}

feature_score_meta#

Meta approach: Meta analysis with each cell deviation contrast result of each feature (a loop for each column in deviation results). And then perform meta analysis with pairwise contrasts to get a pvalue to find to understand what the specific differences are. But there is a problem that meta analysis with so many snps would cost too much time to compute, I am trying to downsample for input, Dan suggest to LD prune with a more permissive threshold.

# Feature score (meta) from per-region contrast results
[feature_score_meta_1]
# File listing per-region posterior-contrast RDS paths
parameter: contrast_units = path
# Random-effects estimator (metafor via pecotmr)
parameter: meta_method = 'REML'
feature_input = [line.split()[0] for line in open(contrast_units).readlines() if line.strip() and not line.strip().startswith('#')]
input: feature_input, group_by = per_chunk
output: f"{cwd}/feature_score_meta/cache/{name}.featurescore{_index+1}.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score.R \
        --method meta \
        --contrast ${_input} \
        --meta-method ${meta_method} \
        --output ${_output}

feature_score_finemap#

Feature score with fine mapped QTL (the top 1 in each CS)

Fine-mapped approach(more recommend by Dan). pick the top SNP in each CS from each cell type fine mapped results. Then pick the most significant one from contrast result. Get the z-score from that snp, which should be the score of that cell type.

# Feature score (finemap) using credible sets from fine-mapping
[feature_score_finemap_1]
parameter: contrast_units = path
# Fine-mapping table RDS with cs_order / pip / variants columns
parameter: fine_mapping = path
# Conditions to score (default: all contexts present in the contrast)
parameter: conditions = ''
feature_input = [line.split()[0] for line in open(contrast_units).readlines() if line.strip() and not line.strip().startswith('#')]
input: feature_input, group_by = per_chunk
output: f"{cwd}/feature_score_finemap/cache/{name}.featurescore{_index+1}.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score.R \
        --method finemap \
        --contrast ${_input} \
        --fine-mapping ${fine_mapping} \
        --conditions "${conditions}" \
        --output ${_output}

# Merge per-chunk finemap feature scores into one table
[feature_score_finemap_2]
input: group_by = "all"
output: f"{cwd}/{name}.{step_name}.feature_score_sum.csv"
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score_merge.R \
        --scores ${_input} \
        --output ${_output}

feature_score_nsig#

Calculate the number of meaningful SNPs in each feature of each cell type (actually, ratio)

# Feature score (nsig) from per-region contrast results
[feature_score_nsig_1]
# File listing per-region posterior-contrast RDS paths
parameter: contrast_units = path
# Significance cutoff for the n-significant ratio
parameter: p_cutoff = 0.00001
feature_input = [line.split()[0] for line in open(contrast_units).readlines() if line.strip() and not line.strip().startswith('#')]
input: feature_input, group_by = per_chunk
output: f"{cwd}/feature_score_nsig/cache/{name}.featurescore{_index+1}.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score.R \
        --method nsig \
        --contrast ${_input} \
        --p-cutoff ${p_cutoff} \
        --output ${_output}

feature_pval_pair#

Find the specific different cell type pair then perform meta analysis with pairwise contrasts to get a pvalue to find to understand what the specific differences are.

# Feature score (pval_pair) from per-region contrast results
[feature_pval_pair_1]
# File listing per-region posterior-contrast RDS paths
parameter: contrast_units = path
parameter: meta_method = 'REML'
# SE floor for the pairwise meta-analysis
parameter: se_cutoff = 0.001
feature_input = [line.split()[0] for line in open(contrast_units).readlines() if line.strip() and not line.strip().startswith('#')]
input: feature_input, group_by = per_chunk
output: f"{cwd}/feature_score_pval_pair/cache/{name}.featurescore{_index+1}.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score.R \
        --method pval_pair \
        --contrast ${_input} \
        --meta-method ${meta_method} \
        --se-cutoff ${se_cutoff} \
        --output ${_output}

# Merge per-chunk feature scores into one table (meta / nsig / pval_pair)
[feature_pval_pair_2, feature_score_meta_2, feature_score_nsig_2]
input: group_by = "all"
output: f"{cwd}/{name}.{step_name}.feature_score_sum.csv"
bash: expand = "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container = container
    Rscript ${modular_script_dir}/pecotmr_integration/mash_feature_score_merge.R \
        --scores ${_input} \
        --output ${_output}