MASH analysis pipeline with data-driven prior matrices

MASH analysis pipeline with data-driven prior matrices#

Fits a MASH model using data-driven prior matrices and computes posterior quantities for the effects of interest.

Overview#

Effects estimated separately in each condition are noisy, and analysing them one condition at a time ignores that most effects are shared. MASH fits a mixture of multivariate normal distributions to the effect estimates, learning from the data which patterns of sharing actually occur - which conditions move together, and how strongly - so that individual estimates can later be shrunk towards those patterns rather than towards zero.

This notebook fits the model. The data-driven prior matrices come from mixture_prior, and applying the fitted model to compute posterior estimates is mash_posterior. By this point the input data have already been converted from the original association summary statistics into the MASH-format object written by mash_preprocessing.

The implementation uses mashr, which improves on the original MASH algorithm in four respects: likelihood and posterior quantities are computed faster through matrix algebra tricks and a C++ implementation; the MASH mixture is fitted faster by convex optimization; the prior can be estimated by approaches other than SFA (see mixture_prior.ipynb); and the estimate of the residual variance \(\hat{V}\) is improved.

The prior matrices come from mixture_prior, fitted in a previous step. After the model is fitted, posteriors are computed for the effects of interest.

When to run it. After mash_preprocessing has assembled effect estimates across conditions and written the MASH-format input file, and before mash_posterior. See mash_preprocessing to generate that file (Urbut et al., 2019).

Input#

--data input/mash/protocol_example.EE.mash.rds (input summary statistics data (MASH-format RDS) produced by mash_preprocessing)

List of 4
 $ random.b: num [1:2000, 1:34] 1.371 -0.565 0.363 0.633 0.404 -0.106 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2000] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ random.s: num [1:2000, 1:34] 1.112 1.056 0.924 1.11 0.996 0.939 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:2000] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ strong.b: num [1:300, 1:34] 1.36 -1.895 1.656 -0.752 1.361 0.64 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ strong.s: num [1:300, 1:34] 1.033 0.966 0.988 0.943 1.212 0.96 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...

--vhat-data input/mash/protocol_example.EE.V_simple.rds (estimated residual correlation matrix across conditions, written by mixture_prior)

 num [1:34, 1:34] 1 0.0625 0.196 0.0857 0.1888 0.0861 ...
 - 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" ...

--prior-data input/mash/protocol_example.EE.prior.rds (data-driven prior covariance matrices U and their weights w, written by mixture_prior)

List of 3
 $ U     :List of 73
  ..$ XtX                         : num [1:34, 1:34] 37.12 -4.7 16.75 -4.48 18.2 -4.56 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ tFLASH_default              : num [1:34, 1:34] 155.052 -0.388 51.932 0.116 56.377 0.443 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_1             : num [1:34, 1:34] 6.06e-03 -7.47e-06 -8.25e-09 -1.81e-05 3.81e-09 -1.02e-05 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_2             : num [1:34, 1:34] 1.06e-02 1.55e-06 3.09e-06 -2.45e-06 8.96e-06 6.54e-07 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_3             : num [1:34, 1:34] 6.06e-03 1.27e-07 3.20e-07 2.81e-07 6.19e-07 2.56e-07 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_4             : num [1:34, 1:34] 6.06e-03 -8.89e-08 -3.04e-08 -1.83e-06 8.94e-09 -2.07e-08 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_5             : num [1:34, 1:34] 5.35e+01 -6.04e-06 1.47e+01 4.66e-01 2.45e+01 4.01e-02 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_6             : num [1:34, 1:34] 6.06e-03 1.16e-09 -8.59e-09 2.99e-08 -1.61e-09 6.04e-10 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_7             : num [1:34, 1:34] 6.06e-03 -5.71e-08 3.58e-08 -7.38e-08 -2.86e-09 2.76e-07 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_8             : num [1:34, 1:34] 0.929 6.251 0.891 2.854 0.858 3.241 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_9             : num [1:34, 1:34] 6.06e-03 1.70e-12 -1.40e-06 1.51e-07 -7.10e-08 1.64e-12 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_10            : num [1:34, 1:34] 6.06e-03 -1.78e-07 3.25e-08 4.60e-08 6.25e-08 -1.30e-06 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_11            : num [1:34, 1:34] 3.69e-02 6.90e-07 -4.80e-02 -3.70e-03 -9.68e-01 1.12e-02 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_12            : num [1:34, 1:34] 5.20e-02 1.35e-04 5.19e-02 4.91e-05 4.40e-02 5.03e-04 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_13            : num [1:34, 1:34] 0.011419 0.000139 0.003632 -0.00019 0.002041 0.000241 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_14            : num [1:34, 1:34] 0.17808 0.00263 0.24694 1.16463 0.09564 2.65106 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_15            : num [1:34, 1:34] 6.06e-03 -4.74e-12 -6.49e-13 1.96e-11 6.18e-11 -1.54e-11 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_16            : num [1:34, 1:34] 0.0413 0.019 0.0209 0.041 0.0178 0.0191 ...
  .. ..- attr(*, "dimnames")=List of 2
  ..$ FLASH_default_17            : num [1:34, 1:34] 6.08e-03 -6.91e-07 6.60e-06 -7.47e-06 -3.12e-06 1.32e-06 ...
  .. ..- attr(*, "dimnames")=List of 2
... (truncated)

--cwd output/mash_fit (working directory all outputs are written under. Defaults to ./mashr_workflow_output.)
--output-prefix protocol_example (prefix for the fitted model and posterior filenames)
--effect-model EE (MASH effect model: EE for exchangeable effects, EZ for exchangeable z-scores)
--compute-posterior True (compute posteriors after fitting. Set False to fit the model only.)

Output#

output/mash_fit/protocol_example_mash.EE.mash_model.rds (fitted MASH model: mixture weights and the fitted prior)

List of 3
 $ mash_model: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
  .. ..$ grid        : num [1:16] 0.0468 0.0663 0.0937 0.1325 0.1874 0.265 ...
  .. ..$ usepointmass: logi TRUE
  ..$ posterior_weights: num [1:2000, 1:18] 0.025 0.0113 0.0268 0.0908 0.015 0.0531 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. ..- attr(*, "names")= chr [1:36000] "1" "604" "650" "700" "723" ...
  ..$ alpha            : num 0
  ..$ lm               :List of 2
  .. ..$ loglik_matrix: num [1:2000, 1:1169] -3.282 -2.661 -6.069 -0.829 -3.036 -1.837 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. ..$ lfactors     : num [1:2000] -45.7 -50.5 -43.9 -43.1 -51.8 -44.3 ...
  ..- attr(*, "class")= chr "mash"
 $ vhat_file : chr <path>/xqtl/jaempawi/xqtl-protocol-new/input/mash/protocol_example.EE.V_simple.rds"
 $ prior_file: chr <path>/xqtl/jaempawi/xqtl-protocol-new/input/mash/protocol_example.EE.prior.rds"

output/mash_fit/protocol_example_mash.EE.posterior.rds (optional; posterior quantities such as lfsr for the “strong” set of gene-SNP pairs)

List of 5
 $ PosteriorMean: num [1:300, 1:34] 0.137 -0.311 0.317 -0.166 0.189 0.117 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ PosteriorSD  : num [1:300, 1:34] 0.344 0.45 0.432 0.47 0.418 0.399 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ lfdr         : num [1:300, 1:34] 0.07816 0.00282 0.03005 0.0025 0.05141 0.03539 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ NegativeProb : num [1:300, 1:34] 0.337 0.744 0.226 0.636 0.313 0.378 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...
 $ lfsr         : num [1:300, 1:34] 0.415 0.256 0.257 0.364 0.364 0.413 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:300] "snp1" "snp2" "snp3" "snp4" "snp5" ...
  .. ..$ : chr [1:34] "ROSMAP_AC_sQTL_PR" "ROSMAP_AC_sQTL_UP" "ROSMAP_DLPFC_sQTL_PR" "ROSMAP_DLPFC_sQTL_UP" "ROSMAP_PCC_sQTL_PR" ...

Minimal Working Example#

Using the MASH prior previously generated, fit the MASH model and compute posteriors.

This example fits the MASH mixture model and computes posteriors for the strong set.

Timing: ~1-2 min (on toy dataset)

sos run pipeline/mash_fit.ipynb mash \
    --output-prefix protocol_example_mash \
    --data input/mash/protocol_example.EE.mash.rds \
    --vhat-data input/mash/protocol_example.EE.V_simple.rds \
    --prior-data input/mash/protocol_example.EE.prior.rds \
    --effect-model EE \
    --compute-posterior \
    --cwd output/mash_fit

Command Interface#

sos run pipeline/mash_fit.ipynb -h

usage: sos run pipeline/mash_fit.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:
  mash

Global Workflow Options:
  --cwd mashr_workflow_output (as path)
  --modular-script-dir code/script (as path)
  --data VAL (as path, required)
                        Input summary statistics data
  --prior-data VAL (as path, required)
  --vhat-data VAL (as path, required)
  --output-prefix ''
                        Prefix of output files. If not specified, it will derive
                        it from data. If it is specified, for example,
                        `--output-prefix AnalysisResults` It will save output
                        files as `{cwd}/AnalysisResults*`.
  --output-suffix all
  --effect-model EE
                        Exchangable effect (EE) or exchangable z-scores (EZ)
  --container ''
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 1 (as int)
                        Number of threads

Sections
  mash_1:               Fit MASH mixture model (learns the mixture weights on
                        the random subset)
    Workflow Options:
      --output-level 4 (as int)
  mash_2:               Compute posterior for the "strong" set of data as in
                        Urbut et al 2017. This is optional because most of the
                        time we want to apply the MASH model learned on much
                        larger data-set.
    Workflow Options:
      --[no-]compute-posterior (default to True)
                        default to True; use --no-compute-posterior to disable
                        this

Workflow implementation#

The mash workflow has two steps. mash_1 fits the mashr mixture model; mash_2 (optional, on by default) computes posteriors for the “strong” set as in Urbut et al 2017. Use --no-compute-posterior to skip the posterior step.

[global]
parameter: cwd = path('./mashr_workflow_output')
parameter: modular_script_dir = path('code/script')
# Input summary statistics data
parameter: data = path
parameter: prior_data = path
parameter: vhat_data = path
# Prefix of output files. If not specified, it will derive it from data.
# If it is specified, for example, `--output-prefix AnalysisResults`
# It will save output files as `{cwd}/AnalysisResults*`.
parameter: output_prefix = ''
parameter: output_suffix = 'all'
# Exchangable effect (EE) or exchangable z-scores (EZ)
parameter: effect_model = 'EE'
parameter: container = ""
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 1
data = data.absolute()
cwd = cwd.absolute()
prior_data = prior_data.absolute()
vhat_data = vhat_data.absolute()
if len(output_prefix) == 0:
    output_prefix = f"{data:bn}"

mash_model = file_target(f"{cwd:a}/{output_prefix}.{effect_model}.mash_model.rds")

`mashr` mixture model fitting#

The main references are the mashr vignettes this for mashr eQTL outline and this for using FLASH prior.

The outcome of this workflow should be found under ./mashr_workflow_output folder (can be configured). File names have pattern *.mash_model_*.rds. They can be used to computer posterior for input list of gene-SNP pairs (see next section).

# Fit MASH mixture model (learns the mixture weights on the random subset)
[mash_1]
parameter: output_level = 4
input: data, vhat_data, prior_data
output: mash_model

task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, 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_fit.R \
        --data ${_input[0]} \
        --vhat-data ${_input[1]} \
        --prior-data ${_input[2]} \
        --effect-model ${effect_model} \
        --output-level ${output_level} \
        --output ${_output}

Optional posterior computations#

Additionally provide posterior for the “strong” set in MASH input data.

# Compute posterior for the "strong" set of data as in Urbut et al 2017.
# This is optional because most of the time we want to apply the
# MASH model learned on much larger data-set.
[mash_2]
# default to True; use --no-compute-posterior to disable this
parameter: compute_posterior = True
# input Vhat file for the batch of posterior data
skip_if(not compute_posterior)

input: data, vhat_data, mash_model
output: f"{cwd:a}/{output_prefix}.{effect_model}.posterior.rds"

task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, 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[0]} \
        --vhat-data ${_input[1]} \
        --mash-model ${_input[2]} \
        --effect-model ${effect_model} \
        --output ${_output}