Help improve this workflow!
This workflow has been published but could be further improved with some additional meta data:- Keyword(s) in categories input, output, operation, topic
You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .
A snakemake workflow to reproduce analysis for the Olink CVD-1 - HF project
:bangbang: Read the full, open-access paper on Circulation
Authors
- Albert Henry (@alhenry)
Usage
The following describes a step-by-step procedure to reproduce the analysis. Code blocks represent input commands for command line interface (e.g. Terminal app on Mac OS)
Step 1: Install Conda & Snakemake
Install Conda / Miniconda (recommended for a fresh installation)
Install Snakemake using:
conda create -c bioconda -c conda-forge -n snakemake snakemake
For installation details, see the instructions in the Snakemake documentation .
Step 2: Obtain a copy of this workflow
Using
git
,
clone
this repository to your local system:
git clone https://github.com/alhenry/cvd1-hf.git
Enter the newly created (cloned) directory:
cd cvd1-hf
Step 3: Execute workflow
Activate the conda environment:
conda activate snakemake
Execute the workflow locally via
snakemake --use-conda -c all
Step 4: Inspect results
A successful execution should create a
results
folder containing 2 files:
res_Observational.csv
(observational results) and
res_MR_aggregate.csv
(MR results)
Notes
Can I run the workflow without snakemake?
Yes, but it requires some tweaks and some familiarity with R.
The main analysis scripts are located in the
workflow/scripts/
folder.
These are standard R scripts embedded in a snakemake pipeline, which can be modified by substituting the relevant snakemake objects with local objects to run specific parts of the analysis.
Change in mr_egger() function in MendelianRandomization package
The sensitivity analysis performed in the present analysis involves running
a MR-Egger regression accounting for correlation between instruments
with
mr_egger(correl = TRUE)
function implemented in
MendelianRandomization
package.
It was recently discovered that there was a minor bug in the
mr_egger(correl = TRUE)
function in which the signs from standardised instrument ratio estimates and correlation matrix are not aligned.
This bug was fixed in the
MendelianRandomization
version
0.4.1
and later, which is now set as the default of this workflow.
For backward reproducibility with earlier version before the fix was introduced,
these settings can be changed by modifying
MR_package_ver
in the
config.yaml
file to an earlier valid version of the package e.g.
0.4.0
(Note that this will force a system-wide installation of the specified version of the package).
Configuring the workflow
As part of the multiverse sensitivity analysis, this workflow is designed to run up to 150 combinations of instrument selection parameters and MR models per protein.
To speed up this computation, the workflow is modularised so that protein-specific analysis can be run separately.
The target proteins for analysis can be configured by modifying the
target_protein_MR
list in the
config.yaml
file (default to all)
Package dependencies
In addition to
conda
,
snakemake
, and
R
, this workflows uses the following
R
packages:
If local installation of R and these packages are already present,
the workflow can be run using local system installation by ommiting the
--use-conda
argument, i.e.:
snakemake -c all
Acknowledgements
Codes for MR with principal component instruments were adapted from Burgess S et al (2017)
Code Snippets
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | source("workflow/scripts/MR_functions.R") # read & reformat data data_MR <- fread(snakemake@input[["data_MR"]]) setnames(data_MR, c("rsID", "beta_Prot", "se_Prot", "beta_HF", "se_HF"), c("SNP", "beta.x", "se.x", "beta.y", "se.y")) snplist <- read_lines(snakemake@input[["snplist"]]) ldrho <- fread(snakemake@input[["ldrho"]], col.names = snplist) %>% as.matrix(rownames.value = snplist) # run MR protein <- snakemake@wildcards[["protein"]] df_res <- data_MR[Protein == protein] %>% group_by(r2_thresh, P_thresh, Protein) %>% nest() %>% # run MR per instrument set mutate(MR = map(data, run_MR_all, ldrho)) %>% select(-data) %>% unnest(cols = MR) fwrite(df_res, snakemake@output[[1]]) |
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | library(tidyverse) library(metafor) # read data df_obs <- read_csv(snakemake@input[[1]]) # helper function to run meta-analysis run_meta <- function(df, y_col = beta, se_col = se, method = "FE", ..., report=c("beta", "se", "pval", "QEp")){ y_col <- enquo(y_col) se_col <- enquo(se_col) betas <- pull(df, !!y_col) ses <- pull(df, !!se_col) meta_res <- tryCatch({ rma(yi = betas, sei = ses, method = method) }, # return NULL if error (e.g. Division by zero) error = function(e) NULL ) df.res <- if (!is.null(meta_res)){ map(report, ~`[[`(meta_res,.) %>% as.numeric) %>% set_names(report) %>% as_tibble } else { rep(NaN, length(report)) %>% as.list %>% set_names(report) %>% as_tibble } df.res } # group dataset per protein df_meta <- df_obs %>% group_by(Protein) %>% nest() %>% # run meta-analysis per-protein mutate(meta = map(data, run_meta)) %>% select(-data) %>% unnest() %>% # convert beta to risk ratio & calculate 95% confidence intervals mutate(RR = exp(beta), RR_LCI = exp(beta - qnorm(0.975)*se), RR_UCI = exp(beta + qnorm(0.975)*se)) # write results write_csv(df_meta, snakemake@output[[1]]) |
20 21 | script: "scripts/Obs_analysis.R" |
32 33 | script: "scripts/MR_analysis.R" |
40 41 | shell: "awk 'NR == FNR {{print; next}} FNR > 1 {{print}}' {input} > {output}" |
Support
- Share:
-
-
-
- Future updates
Related Workflows
