In jwr-git/mrpipeline: Implements a Ready-for-Use Mendelian Randomisation Pipeline

knitr::opts_chunk$set(echo = F)

Report Overview

This page contains an overview for the analyses that have been carried out. Links are provided to each per-trait specific analysis in the dataset tables. All tables and figures are downloadable and all separate report pages are self-contained html files, meaning these can be sent to other people without extra files.

Details and explanations of analyses and results are given throughout. Certain terms related to MR are hyperlinked to the MR Dictionary, which offers an explanation for that term.

Overview of Results

This shows a basic overview of all of the results, which are given in more detail below and in individual reports. This overview shows a basic "Pass" or "Fail" mark if that exposure-outcome pair passes or fails that analysis. Passes are defined for each of the analyses:

• MR: P value < r conf$mr_p_thres.
• Steiger filtering: P < 0.05 and correctly estimated direction of effect.
• Colocalisation: H4 > r conf$coloc_strong_thres is strong evidence, H4 > r conf$coloc_moderate_thres is moderate evidence, and weak evidence otherwise.

mrres <- params$report$results %>%
  dplyr::select(id.exposure, id.outcome, exposure, outcome, pval) %>%
  dplyr::mutate(mrflag = ifelse(pval < conf$mr_p_thres, "Pass", "Fail"))

steigerres <- params$report$steigerresults %>%
  dplyr::select(id.exposure, id.outcome, flag) %>%
  dplyr::mutate(flag = ifelse(flag == "True", "Pass",
                              ifelse(flag == "Unknown", "Uncertain", "Fail")))

colocres <- params$report$cresults %>%
  dplyr::select(name1, name2, nsnps, H4) %>%
  dplyr::rename(exposure = name1, outcome = name2) %>%
  dplyr::mutate(colocflag = ifelse(H4 > conf$coloc_strong_thres & nsnps > 50, "Strong",
                                   ifelse(H4 > conf$coloc_moderate_thres & nsnps > 50, "Moderate", "Weak"))) %>%
  dplyr::mutate(colocflag = ifelse(colocflag == "", "NA", colocflag))

dgidbres <- params$report$dgidbresults %>%
  dplyr::select(id, Druggable.Genome, Clinically.Actionable, Drug.Resistant) %>%
  dplyr::mutate(Ontologies = paste0(
    ifelse(Druggable.Genome == T, "Drg.Gen", ""),
    ifelse(Clinically.Actionable == T, "|Clin.Act", ""),
    ifelse(Drug.Resistant == T, "|Drg.Res", "")
  ))

drugs <- params$report$dgidbdrugs %>%
  dplyr::count(id)

dplyr::left_join(mrres, steigerres, by = c("id.exposure", "id.outcome")) %>%
  dplyr::left_join(colocres, by = c("exposure", "outcome")) %>%
  dplyr::left_join(dgidbres[c("id", "Ontologies")], by = c("id.exposure" = "id")) %>%
  dplyr::left_join(drugs, by = c("id.exposure" = "id")) %>%
  dplyr::select(-pval, -nsnps, -H4) %>%
  dplyr::rename(Exposure.ID = id.exposure,
                Outcome.ID = id.outcome,
                Exposure = exposure,
                Outcome = outcome,
                MR.Pass = mrflag,
                Steiger.Pass = flag,
                Coloc.Pass = colocflag,
                Drugs.Count = n) %>%
  dplyr::mutate(Exposure = paste0("<a href='./traits/", Exposure, ".html'>", Exposure, "</a>")) %>%
  dplyr::mutate(Outcome = paste0("<a href='./traits/", Outcome, ".html'>", Outcome, "</a>")) %>%
  dplyr::mutate(dplyr::across(Drugs.Count, ~tidyr::replace_na(., 0))) %>%
  dplyr::relocate(Exposure, Outcome,
                  Exposure.ID, Outcome.ID,
                  MR.Pass, Steiger.Pass, Coloc.Pass) %>%
  DT::datatable(rownames = F,
                escape = F,
                caption = "Overview of all results",
                filter = "top",
                extensions = c("Buttons"),
                options = list(dom = "BDfrtip", 
                               buttons = list('copy',
                                              list(extend = 'collection', 
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download')))) %>%
  DT::formatStyle(
    "MR.Pass",
    color = DT::styleEqual(c("Pass", "Fail"),
                           c("darkgreen", "red")),
    fontWeight = DT::styleEqual("Pass", "bold")
  ) %>%
  DT::formatStyle(
    "Steiger.Pass",
    color = DT::styleEqual(c("Pass", "Uncertain", "Fail"),
                           c("darkgreen", "orange", "red")),
    fontWeight = DT::styleEqual("Pass", "bold")
  ) %>%
  DT::formatStyle(
    "Coloc.Pass",
    color = DT::styleEqual(c("Strong", "Moderate", "Weak", "NA"),
                           c("darkgreen", "orange", "red", "red")),
    fontWeight = DT::styleEqual("Strong", "bold")
  )# %>% 
#  DT::formatStyle(
#    c("Drug.Genome", "Clin.Act", "Drug.Resistant"),
#    color = DT::styleEqual(c("Y", "N"),
#                           c("darkgreen", "red")),
#    fontWeight = DT::styleEqual("Y", "bold")
#  )

Datasets

Tables presented here contain a breakdown and overview of the datasets. Trait.Name links take you to the specific trait results page and the Trait.ID links will take you to the OpenGWAS DB report, if applicable for that trait.

Exposure datasets

There were r nrow(unique(params$report$exposures)) exposure datasets provided. Details are provided in the following table for each dataset, including how many SNPs were removed at each stage of the QC pipeline.

params$report$exposures %>%
  dplyr::mutate(Trait.Annotated = paste0("<a href='./traits/", Trait.Annotated, ".html'>", Trait.Annotated, "</a>"),
                Trait.ID = ifelse(From.IEUGWASDB == 1, paste0("<a href='https://gwas.mrcieu.ac.uk/datasets/", Trait.ID, "/'>", Trait.ID, "</a>"), Trait.ID)) %>%
  dplyr::relocate(Trait.ID, .after = Trait.Annotated) %>%
  dplyr::select(-From.IEUGWASDB, -Trait.Name) %>%
  DT::datatable(rownames = F,
                escape = F,
                filter = 'top',
                extensions = 'Buttons',
                caption = "Breakdown for each exposure dataset.",
                options = list(dom = 'Bfrtip',
                               buttons = list('copy',
                                              list(extend = 'collection',
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download'))))

Column headers and meanings

• Trait.Annotated: Name of the trait - Links to that trait's results report.
• Trait.ID: Trait ID - Links to the OpenGWAS DB, if applicable.
• SNPs: Number of SNPs initially loaded from the GWAS (either from the file or from the OpenGWAS database).
• Preclumped: Whether or not the dataset was already clumped (e.g., data from OpenGWAS can be retrieved after undergoing clumping with the standard parameters).
• Rem.F.Stat: Number of SNPs which were removed due to not meeting the F-statistic cut-off.
• Rem.Clumped: Number of SNPs which were removed due to clumping. Note that this will be 0 if the data was pre-clumped.
• SNPs.Formatted: Number of SNPs which were successfully formatted to the MR-Base/TwoSampleMR package format (e.g., SNPs which had negative standard errors are removed at this stage).

Exposure SNPs

This table details the final list of formatted SNPs for the exposure datasets which were used in the MR analyses.

params$dat[, colnames(params$dat) %in% c("SNP",
                                        "effect_allele.exposure",
                                        "other_allele.exposure",
                                        "beta.exposure",
                                        "eaf.exposure",
                                        "chr",
                                        "pos",
                                        "se.exposure",
                                        "pval.exposure",
                                        "exposure",
                                        "mr_keep.exposure",
                                        "f.stat.exposure")] %>%
  unique() %>%
  dplyr::rename(Eff.Allele = effect_allele.exposure,
                Alt.Allele = other_allele.exposure,
                Beta = beta.exposure,
                EAF = eaf.exposure,
                Chr = chr,
                Pos = pos,
                SE = se.exposure,
                P.value = pval.exposure,
                Exposure = exposure,
                Inc.MR = mr_keep.exposure,
                F.Stat = f.stat.exposure) %>%
  DT::datatable(rownames = F,
                caption = "Details for SNPs which were included in the MR analysis by passing the QC and cleaning stages.",
                extensions = c("Buttons"),
                options = list(order = list(8, "desc"), 
                               dom = "BDfrtip",
                               buttons = list('copy',
                                              list(extend = 'collection', 
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download')))) %>%
  DT::formatSignif(columns = c("P.value"), digits = 3) %>%
  DT::formatRound(columns = c("F.Stat"), digits = 0)

Outcome datasets

There were r nrow(unique(params$report$outcomes)) outcome datasets provided. Details are provided in the following table for each dataset, including details on how many SNPs were proxied.

params$report$outcomes %>%
  dplyr::mutate(Trait.Annotated = paste0("<a href='./traits/", Trait.Annotated, ".html'>", Trait.Annotated, "</a>"),
                Trait.ID = ifelse(From.IEUGWASDB == 1, paste0("<a href='https://gwas.mrcieu.ac.uk/datasets/", Trait.ID, "/'>", Trait.ID, "</a>"), Trait.ID)) %>%
  dplyr::relocate(Trait.ID, .after = Trait.Annotated) %>%
  dplyr::select(-From.IEUGWASDB, -Trait.Name) %>%
  DT::datatable(rownames = F,
                escape = F,
                filter = 'top',
                extensions = 'Buttons',
                caption = "Breakdown for each outcome dataset.",
                options = list(dom = 'Bfrtip',                                
                               buttons = list('copy',
                                              list(extend = 'collection', 
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download'))))

Column headers and meanings

• Trait.Annotated: Name of the trait - Links to that trait's results report.
• Trait.ID: Trait ID - Links to the OpenGWAS DB, if applicable.
• SNPs: Number of SNPs initially loaded from the GWAS (either from the file or from the OpenGWAS database).
• Num.Proxies: Number of proxied SNPs. A list of all unique SNPs from all exposure datasets are used to extract outcome SNPs -- if a SNP is not found then a proxy SNP is used instead.
• SNPs.Formatted: Number of SNPs which were successfully formatted to the MR-Base/TwoSampleMR package format (e.g., SNPs which had negative standard errors are removed at this stage).

Outcome SNPs

This table details the final list of formatted SNPs for the outcome datasets which were used in the MR analyses.

params$dat[, colnames(params$dat) %in% c("SNP",
                                        "effect_allele.outcome",
                                        "other_allele.outcome",
                                        "beta.outcome",
                                        "eaf.outcome",
                                        "chr",
                                        "pos",
                                        "se.outcome",
                                        "pval.outcome",
                                        "outcome",
                                        "proxy.outcome",
                                        "proxy_snp.outcome")] %>%
  dplyr::rename(Eff.Allele = effect_allele.outcome,
                Alt.Allele = other_allele.outcome,
                Beta = beta.outcome,
                EAF = eaf.outcome,
                Chr = chr,
                Pos = pos,
                SE = se.outcome,
                P.value = pval.outcome,
                Outcome = outcome,
                Proxy = proxy.outcome,
                Proxy.SNP = proxy_snp.outcome) %>%
  DT::datatable(rownames = F,
                caption = "Details for outcome SNPs, or their proxies, which were included in the MR analysis.",
                extensions = c("Buttons"),
                options = list(order = list(8, "asc"), 
                               dom = "BDfrtip", 
                               buttons = list('copy',
                                              list(extend = 'collection', 
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download')))) %>%
  DT::formatSignif(columns = c("P.value"), digits = 3)

MR Results

MR results are presented here for the Wald ratio (WR) and inverse variance weighted (IVW). The per-trait reports detail further MR methods, such as MR-Egger, but are omitted here for brevity.

There were a total of r 0.05 / params$report$bonferroni MR tests conducted. A Bonferroni P-value threshold would therefore be r params$report$bonferroni ($0.05 / tests$).

steiger <- params$report$steigerresults %>% 
  dplyr::select(id.exposure, id.outcome, flag)

params$report$results %>%
  dplyr::left_join(steiger) %>%
  dplyr::rename(Exposure.ID = id.exposure,
                Outcome.ID = id.outcome,
                Outcome = outcome,
                Exposure = exposure,
                Method = method,
                SNPs = nsnp,
                P.value = pval,
                OR = or,
                Lower.95ci = or_lci95,
                Upper.95ci = or_uci95,
                Steiger.Flag = flag) %>%
  dplyr::select(-Exposure.ID, -Outcome.ID) %>%
  dplyr::relocate(Outcome, .after = Exposure) %>%
  DT::datatable(rownames = F,
                filter = 'top',
                extensions = c("Buttons"),
                options = list(order = list(4, "asc"), 
                               dom = "Bfrtip", 
                               buttons = list('copy',
                                              list(extend = 'collection', 
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download')))) %>%
  DT::formatSignif(columns = c("P.value", "OR", "Lower.95ci", "Upper.95ci"), digits = 3)

Column headers and meanings

• Exposure: Exposure trait names.
• Outcome: Outcome trait names.
• Method: MR method, either WR for single-SNP instruments or IVW for multi-SNP instruments.
• SNPs: Number of SNPs that form the instrument.
• P.value: MR P value.
• OR: MR odds ratios.
• Lower.95ci: Lower 95% confidence interval.
• Upper.95ci: Upper 95% confidence interval.
• Steiger.Flag: Flag for the Steiger filtering results: true, direction of effect is from exposure to the outcome; false, the direction of effect is from outcome to the exposure; unknown, Steiger P value is > 0.05.

Colocalisation Results

Here are all of the colocalisation results. Regional plots are available in the per-trait reports.

if (length(params$report$plots[params$report$plots$type == "regional", ]) > 0) {
  for (idx in as.integer(rownames(params$report$plots[params$report$plots$type == "regional", ])))
  {
    if (all(is.na(params$report$raw.plots[[idx]]))) {
      next
    }
    grid.newpage()
    grid.draw(params$report$raw.plots[[idx]])
  }
}

cresults <- params$report$cresults

for (i in 1:nrow(cresults)) {
    tryCatch({
      cresults[i, "plot"] <- sprintf("<a class='ItemsTooltip' target='_blank'><img class='imgTooltip' src='%s'/>Region</a>", knitr::image_uri(knitr::fig_chunk('regional-plot', 'png', i)))
    },
    error = function(e) {
      cresults[i, "plot"] <- "No plot"
    })
}
cresults <- cresults[cresults["nsnps"] > 50, ]

if (length(cresults)) {
  cresults %>%
    dplyr::rename(Exposure = name1,
                  Outcome = name2,
                  Plot = plot,
                  SNPs = nsnps,
                  ChrPos = chrpos) %>%
    dplyr::select(-id1, -id2) %>%
    dplyr::relocate(Plot, .after = Outcome) %>%
    DT::datatable(rownames = F,
                  escape = F,
                  filter = 'top',
                  extensions = c("Buttons"),
                  options = list(order = list(8, "desc"), 
                                 dom = "Bfrtip", 
                                 buttons = list('copy',
                                              list(extend = 'collection', 
                                                   buttons = c('csv', 'excel'), 
                                                   text = 'Download')))) %>%
    DT::formatPercentage(c("H0", "H1", "H2", "H3", "H4"), 2)
} else {
  cat("There are no colocalisation results.")
}

Column headers and meanings

• Exposure: Name for Exposures.
• Outcome: Name for Outcomes.
• Plot: Hover over for the regional plot.
• SNPs: Number of SNPs that form the instrument.
• H0-4: Posterior probability for hypotheses 0 through 4.
• ChrPos: Chromosome and position range used to extract SNPs for colocalisation.

Overview of Drug Target Evidence

# brks <- seq(from = 0, to = 1, by = 0.1)
# clrs <- round(seq(180, 40, length.out = length(brks) + 1), 0) %>%
#   {paste0("rgb(", ., ",", ., ",180)")}
# 
# to_table <- params$report$otresults %>%
#   dplyr::select(-id1, -id2) %>%
#   dplyr::relocate(name2, .after = name1) %>%
#   dplyr::rename(Exposure = name1,
#                 Outcome = name2,
#                 Overall.Score = overall.ot,
#                 Literature = literature.ot,
#                 RNA.Expression = rna_expression.ot,
#                 Genetic.Assoc = genetic_assoc.ot,
#                 Somatic.Mutations = somatic_mute.ot,
#                 Known.Drugs = known_drug.ot,
#                 Animal.Model = animal_model.ot,
#                 Affected.Pathways = affected_pathway.ot)
# 
# to_table %>%
#   DT::datatable(rownames = F,
#                 escape = F,
#                 filter = 'top',
#                 extensions = c("Buttons"),
#                 options = list(order = list(2, "desc"), 
#                                dom = "Bfrtip", 
#                                buttons = list('copy',
#                                               list(extend = 'collection', 
#                                                    buttons = c('csv', 'excel'), 
#                                                    text = 'Download')))) %>%
#     DT::formatSignif(columns = c("Overall.Score", 
#                                "Literature", 
#                                "RNA.Expression", 
#                                "Genetic.Assoc",
#                                "Somatic.Mutations",
#                                "Known.Drugs",
#                                "Animal.Model",
#                                "Affected.Pathways"), 
#                    digits = 3) %>%
#   DT::formatStyle(names(to_table),
#                   background = DT::styleColorBar(range(0, 1), 'lightblue'),
#                   backgroundSize = '98% 88%',
#                   backgroundRepeat = 'no-repeat',
#                   backgroundPosition = 'center')

Literature Evidence

Evidence from the literature collected by EpigraphDB.

if (!any(is.null(params$report$epidbresults)) & nrow(params$report$epidbresults) > 0) {
  to_table <- params$report$epidbresults
  to_table$pubmed_id <- lapply(to_table$pubmed_id, function(x) paste(lapply(unlist(x), function(y) pmid_to_link(y, hyperlink = T)), collapse = ", "))
  #to_table$pubmed_id <- gsub(",", ", ", to_table$pubmed_id)
  to_table %>%
    dplyr::select(-id1, -id2, -lt.id, -lt.type, -literature_count) %>%
    dplyr::relocate(c("gene.name", "st.predicate", "lt.name", "pubmed_id")) %>%
    dplyr::rename(Gene = gene.name,
                  PMIDs = pubmed_id,
                  Outcome = lt.name,
                  Predicate = st.predicate) %>%
    DT::datatable(rownames = F,
                  escape = F,
                  filter = 'top',
                  extensions = c("Buttons"),
                  options = list(order = list(0, "asc"),
                                 dom = "Bfrtip",
                                 buttons = list ('copy',
                                                 list(extend = 'collection',
                                                      buttons = c('csv', 'excel'),
                                                      text = 'Download'))))
}

r if (all(is.na(params$report$dgidbresults))) { "\\begin{comment}" }

Drug-Genome Interaction DB

params$report$dgidbresults %>%
  dplyr::select(-id) %>%
  dplyr::rename(Trait = trait) %>%
  dplyr::mutate(Druggable.Genome = ifelse(Druggable.Genome == T, "Yes", "No"),
                Clinically.Actionable = ifelse(Clinically.Actionable == T, "Yes", "No"),
                Drug.Resistant = ifelse(Drug.Resistant == T, "Yes", "No")) %>%
  DT::datatable(rownames = F,
                escape = F,
                filter = 'top',
                extensions = c("Buttons"),
                options = list(order = list(0, "asc"),
                               dom = "Bfrtip",
                               buttons = list ('copy',
                                               list(extend = 'collection',
                                                    buttons = c('csv', 'excel'),
                                                    text = 'Download')))) %>%
  DT::formatStyle(
    c("Druggable.Genome", "Clinically.Actionable", "Drug.Resistant"),
    color = DT::styleEqual(c("Yes", "No"),
                           c("darkgreen", "red")),
    fontWeight = DT::styleEqual("Yes", "bold")
  )

Possible Drugs

to_table <- params$report$dgidbdrugs
to_table$PMIDs <- lapply(to_table$PMIDs, function(x) paste(lapply(unlist(x), function(y) pmid_to_link(y, hyperlink = T)), collapse = ", "))

to_table %>%
  dplyr::select(-id) %>%
  dplyr::rename(Trait = trait) %>%
  DT::datatable(rownames = F,
                escape = T,
                filter = 'top',
                extensions = c("Buttons"),
                options = list(dom = "Bfrtip",
                               buttons = list ('copy',
                                               list(extend = 'collection',
                                                    buttons = c('csv', 'excel'),
                                                    text = 'Download'))))

r if (all(is.na(params$report$dgidbresults))) { "\\end{comment}" }

jwr-git/mrpipeline documentation built on Oct. 2, 2022, 3:41 p.m.