This RMarkdown document demonstrates how key elements from the notebook for case study 1 in the EpiGraphDB paper can be achieved using the R package. For detailed explanations of the case study please refer to the paper or the case study notebook.

Context

Mendelian randomization (MR), a technique to evaluate the causal role of modifiable exposures on a health outcome using a set of genetic instruments, tacitly assumes that a genetic variant (e.g. SNP) is only related to an outcome of interest through an exposure (i.e. the “exclusion restriction criterion”). Horizontal pleiotropy, where a SNP is associated with multiple phenotypes independently of the exposure of interest, potentially violates this assumption delivering misleading conclusions. In contrast, vertical pleiotropy, where a SNP is associated with multiple phenotypes on the same biological pathway, does not violate this assumption.

Here, we use external evidence on biological pathways and protein-protein interactions to assess the pleiotropic profile for a genetic variant based on the genes (and proteins) with which that variant is associated. In a graph representation, we will show that:

  • We have evidence of potential vertical pleiotropy when the proteins associated with a variant belong to one connected group,
  • In contrast, if the proteins associated with a single variant are isolated from each other, this suggests evidence of horizontal pleiotropy.

This case study goes as follows:

  • For a list of genes that we wish to assess their pleiotropic profiles (against a variant of interests), we retrieve from EpiGraphDB their mapped protein products.
  • Then respectively from the evidence on shared pathway network and protein-protein interaction network, we can extract the proteins groups formed in these networks as connected communities to assess the pleiotropic profiles of the genes.

Here we configure the parameters used in the case study example. In this example we will look at the variant rs12720356, a SNP located in chromosome 19 that has been associated with Crohn’s disease and psoriasis.**

SNP <- "rs12720356"

GENELIST <- c("ZGLP1", "FDX1L", "MRPL4", "ICAM5", "TYK2", "GRAP2", "KRI1", "TMED1", "ICAM1")

PPI_N_INTERMEDIATE_PROTEINS <- 1

Mapping genes to proteins

The first step of the analysis is to map each of these genes to their protein product.

get_gene_protein <- function(genelist) {
  endpoint <- "/mappings/gene-to-protein"
  params <- list(
    gene_name_list = genelist %>% I()
  )
  r <- query_epigraphdb(
    route = endpoint,
    params = params,
    mode = "table",
    method = "POST"
  )
  protein_df <- r
  if (nrow(protein_df) > 0) {
    res_df <- protein_df %>%
      select(gene_name = `gene.name`, uniprot_id = `protein.uniprot_id`)
  } else {
    res_df <- tibble() %>% set_names(c("gene_name", "uniprot_id"))
  }

  res_df
}

gene_protein_df <- get_gene_protein(genelist = GENELIST)
gene_protein_df
#> # A tibble: 9 x 2
#>   gene_name uniprot_id
#>   <chr>     <chr>     
#> 1 ZGLP1     P0C6A0    
#> 2 FDX1L     Q6P4F2    
#> 3 MRPL4     Q9BYD3    
#> 4 ICAM5     Q9UMF0    
#> 5 TYK2      P29597    
#> 6 GRAP2     O75791    
#> 7 KRI1      Q8N9T8    
#> 8 TMED1     Q13445    
#> 9 ICAM1     P05362

Identifying the involvement in the same biological pathways

For each protein we retrieve the pathways they are involved in.

get_protein_pathway <- function(gene_protein_df) {
  endpoint <- "/protein/in-pathway"
  params <- list(
    uniprot_id_list = gene_protein_df %>% pull(`uniprot_id`) %>% I()
  )
  df <- query_epigraphdb(route = endpoint, params = params, mode = "table", method = "POST")

  if (nrow(df) > 0) {
    res_df <- gene_protein_df %>%
      select(`uniprot_id`) %>%
      left_join(df, by = c("uniprot_id"))
  } else {
    res_df <- gene_protein_df %>%
      select(`uniprot_id`) %>%
      mutate(pathway_count = NA_integer_, pathway_reactome_id = NA_character_)
  }
  res_df <- res_df %>%
    mutate(
      pathway_count = ifelse(is.na(pathway_count), 0L, as.integer(pathway_count)),
      pathway_reactome_id = ifelse(is.na(pathway_reactome_id), c(), pathway_reactome_id)
    )

  res_df
}
pathway_df <- get_protein_pathway(gene_protein_df = gene_protein_df)
pathway_df
#> # A tibble: 9 x 3
#>   uniprot_id pathway_count pathway_reactome_id
#>   <chr>              <int> <list>             
#> 1 P0C6A0                 0 <NULL>             
#> 2 Q6P4F2                15 <chr [15]>         
#> 3 Q9BYD3                 6 <chr [6]>          
#> 4 Q9UMF0                 5 <chr [5]>          
#> 5 P29597                28 <chr [28]>         
#> 6 O75791                15 <chr [15]>         
#> 7 Q8N9T8                 0 <NULL>             
#> 8 Q13445                 0 <NULL>             
#> 9 P05362                11 <chr [11]>

Now for each pair of proteins we match the pathways they have in common.

get_shared_pathway <- function(pathway_df) {
  # For the protein-pathway data
  # Get protein-protein permutations where they share pathways

  per_permutation <- function(pathway_df, permutation) {
    df <- pathway_df %>% filter(uniprot_id %in% permutation)
    primary_pathway <- pathway_df %>%
      filter(uniprot_id == permutation[1]) %>%
      pull(pathway_reactome_id) %>%
      unlist()
    assoc_pathway <- pathway_df %>%
      filter(uniprot_id == permutation[2]) %>%
      pull(pathway_reactome_id) %>%
      unlist()
    intersect(primary_pathway, assoc_pathway)
  }

  pairwise_permutations <- pathway_df %>%
    pull(`uniprot_id`) %>%
    gtools::permutations(n = length(.), r = 2, v = .)
  shared_pathway_df <- tibble(
    protein = pairwise_permutations[, 1],
    assoc_protein = pairwise_permutations[, 2]
  ) %>%
    mutate(
      shared_pathway = map2(`protein`, `assoc_protein`, function(x, y) {
        per_permutation(pathway_df = pathway_df, permutation = c(x, y))
      }),
      combination = map2_chr(`protein`, `assoc_protein`, function(x, y) {
        comb <- sort(c(x, y))
        paste(comb, collapse = ",")
      }),
      count = map_int(`shared_pathway`, function(x) na.omit(x) %>% length()),
      connected = count > 0
    )

  shared_pathway_df
}
shared_pathway_df <- get_shared_pathway(pathway_df)
n_pairs <- length(shared_pathway_df %>% filter(count > 0))
print(glue::glue("Num. shared_pathway pairs: {n_pairs}"))
#> Num. shared_pathway pairs: 6
shared_pathway_df %>% arrange(desc(count))
#> # A tibble: 72 x 6
#>    protein assoc_protein shared_pathway combination   count connected
#>    <chr>   <chr>         <list>         <chr>         <int> <lgl>    
#>  1 P05362  P29597        <chr [6]>      P05362,P29597     6 TRUE     
#>  2 P29597  P05362        <chr [6]>      P05362,P29597     6 TRUE     
#>  3 P05362  Q9UMF0        <chr [5]>      P05362,Q9UMF0     5 TRUE     
#>  4 Q9UMF0  P05362        <chr [5]>      P05362,Q9UMF0     5 TRUE     
#>  5 O75791  P05362        <chr [2]>      O75791,P05362     2 TRUE     
#>  6 O75791  P29597        <chr [2]>      O75791,P29597     2 TRUE     
#>  7 O75791  Q9UMF0        <chr [2]>      O75791,Q9UMF0     2 TRUE     
#>  8 P05362  O75791        <chr [2]>      O75791,P05362     2 TRUE     
#>  9 P29597  O75791        <chr [2]>      O75791,P29597     2 TRUE     
#> 10 Q9UMF0  O75791        <chr [2]>      O75791,Q9UMF0     2 TRUE     
#> # … with 62 more rows

We can further query EpiGraphDB regarding the detailed pathway information using GET /meta/nodes/Pathway/search.

get_pathway_info <- function(reactome_id) {
  endpoint <- "/meta/nodes/Pathway/search"
  params <- list(id = reactome_id)
  df <- query_epigraphdb(route = endpoint, params = params, mode = "table")
  df
}

pathway <- shared_pathway_df %>%
  pull(shared_pathway) %>%
  unlist() %>%
  unique()

pathway_info <- pathway %>% map_df(get_pathway_info)
pathway_info %>% print()
#> # A tibble: 12 x 6
#>    node._name        node.name        node._source node.id node._id node.url    
#>    <chr>             <chr>            <list>       <chr>   <chr>    <chr>       
#>  1 Immune System     Immune System    <chr [1]>    R-HSA-… R-HSA-1… https://rea…
#>  2 Adaptive Immune … Adaptive Immune… <chr [1]>    R-HSA-… R-HSA-1… https://rea…
#>  3 Signal Transduct… Signal Transduc… <chr [1]>    R-HSA-… R-HSA-1… https://rea…
#>  4 Interferon Signa… Interferon Sign… <chr [1]>    R-HSA-… R-HSA-9… https://rea…
#>  5 Interleukin-4 an… Interleukin-4 a… <chr [1]>    R-HSA-… R-HSA-6… https://rea…
#>  6 Interleukin-10 s… Interleukin-10 … <chr [1]>    R-HSA-… R-HSA-6… https://rea…
#>  7 Signaling by Int… Signaling by In… <chr [1]>    R-HSA-… R-HSA-4… https://rea…
#>  8 Cytokine Signali… Cytokine Signal… <chr [1]>    R-HSA-… R-HSA-1… https://rea…
#>  9 Integrin cell su… Integrin cell s… <chr [1]>    R-HSA-… R-HSA-2… https://rea…
#> 10 Immunoregulatory… Immunoregulator… <chr [1]>    R-HSA-… R-HSA-1… https://rea…
#> 11 Extracellular ma… Extracellular m… <chr [1]>    R-HSA-… R-HSA-1… https://rea…
#> 12 Disease           Disease          <chr [1]>    R-HSA-… R-HSA-1… https://rea…

In order to extract protein groups from the shared pathways, the last step for this query is to convert the shared pathway data into a graph where

  • Proteins are represented as nodes.
  • And shared pathways are represented as edges.
  • Then the protein groups (via shared pathways) are the connected communities in the graph.

We then count the number of nodes in each connected community and plot the graph.

protein_df_to_graph <- function(df) {
  df_connected <- df %>%
    filter(connected) %>%
    distinct(`combination`, .keep_all = TRUE)
  nodes <- df %>%
    pull(protein) %>%
    unique()
  graph <- igraph::graph_from_data_frame(
    df_connected,
    directed = FALSE, vertices = nodes
  )
  graph$layout <- igraph::layout_with_kk
  graph
}

graph_to_protein_groups <- function(graph) {
  graph %>%
    igraph::components() %>%
    igraph::groups() %>%
    tibble(group_member = .) %>%
    mutate(group_size = map_int(`group_member`, length)) %>%
    arrange(desc(group_size))
}

pathway_protein_graph <- shared_pathway_df %>% protein_df_to_graph()
pathway_protein_groups <- pathway_protein_graph %>% graph_to_protein_groups()
pathway_protein_groups %>% str()
#> tibble [5 × 2] (S3: tbl_df/tbl/data.frame)
#>  $ group_member:List of 5
#>   ..$ 1: chr [1:5] "O75791" "P05362" "P29597" "Q6P4F2" ...
#>   ..$ 2: chr "P0C6A0"
#>   ..$ 3: chr "Q13445"
#>   ..$ 4: chr "Q8N9T8"
#>   ..$ 5: chr "Q9BYD3"
#>   ..- attr(*, "dim")= int 5
#>   ..- attr(*, "dimnames")=List of 1
#>   .. ..$ : chr [1:5] "1" "2" "3" "4" ...
#>  $ group_size  : Named int [1:5] 5 1 1 1 1
#>   ..- attr(*, "names")= chr [1:5] "1" "2" "3" "4" ...
plot(pathway_protein_graph)

Retrieving shared protein-protein interactions

We can complement the information given by pathway ontologies with data of interactions between proteins. For each pair of proteins, we retrieve shared protein-protein interactions, either a direct interaction or an interaction where there’s one mediator protein in the middle (see parameters section).

get_ppi <- function(gene_protein_df, n_intermediate_proteins = 0) {
  endpoint <- "/protein/ppi/pairwise"
  params <- list(
    uniprot_id_list = gene_protein_df %>% pull("uniprot_id") %>% I(),
    n_intermediate_proteins = n_intermediate_proteins
  )
  df <- query_epigraphdb(
    route = endpoint,
    params = params,
    mode = "table",
    method = "POST"
  )

  if (nrow(df) > 0) {
    res_df <- gene_protein_df %>%
      select(protein = uniprot_id) %>%
      left_join(df, by = c("protein"))
  } else {
    res_df <- gene_protein_df %>%
      select(protein = uniprot_id) %>%
      mutate(assoc_protein = NA_character_, path_size = NA_integer_)
  }
  res_df
}
ppi_df <- get_ppi(gene_protein_df, n_intermediate_proteins = PPI_N_INTERMEDIATE_PROTEINS)
ppi_df
#> # A tibble: 1,686 x 3
#>    protein assoc_protein path_size
#>    <chr>   <chr>             <int>
#>  1 P0C6A0  P0C6A0                2
#>  2 P0C6A0  P0C6A0                2
#>  3 Q6P4F2  Q6P4F2                2
#>  4 Q6P4F2  Q6P4F2                2
#>  5 Q6P4F2  Q6P4F2                2
#>  6 Q6P4F2  Q6P4F2                2
#>  7 Q6P4F2  Q6P4F2                2
#>  8 Q6P4F2  Q6P4F2                2
#>  9 Q6P4F2  Q6P4F2                2
#> 10 Q6P4F2  Q6P4F2                2
#> # … with 1,676 more rows

Like in the query for pathway data, we convert and analyse these results as a graph.

get_shared_ppi <- function(ppi_df) {
  per_row <- function(ppi_df, x, y) {
    res <- ppi_df %>%
      filter(
        `protein` == x,
        `assoc_protein` == y
      ) %>%
      nrow() %>%
      `>`(0)
    res
  }

  pairwise_permutations <- ppi_df %>%
    pull(`protein`) %>%
    unique() %>%
    gtools::permutations(n = length(.), r = 2, v = .)
  shared_ppi_df <- tibble(
    protein = pairwise_permutations[, 1],
    assoc_protein = pairwise_permutations[, 2]
  ) %>%
    mutate(
      shared_ppi = map2_lgl(`protein`, `assoc_protein`, function(x, y) {
        per_row(ppi_df = ppi_df, x = x, y = y)
      }),
      combination = map2_chr(`protein`, `assoc_protein`, function(x, y) {
        comb <- sort(c(x, y))
        paste(comb, collapse = ",")
      }),
      connected = shared_ppi
    )

  shared_ppi_df
}

shared_ppi_df <- get_shared_ppi(ppi_df)
n_ppi_pairs <- shared_ppi_df %>%
  filter(shared_ppi) %>%
  nrow()
print(glue::glue("Num. shared_ppi pairs: {n_ppi_pairs}"))
#> Num. shared_ppi pairs: 10
shared_ppi_df
#> # A tibble: 72 x 5
#>    protein assoc_protein shared_ppi combination   connected
#>    <chr>   <chr>         <lgl>      <chr>         <lgl>    
#>  1 O75791  P05362        TRUE       O75791,P05362 TRUE     
#>  2 O75791  P0C6A0        FALSE      O75791,P0C6A0 FALSE    
#>  3 O75791  P29597        TRUE       O75791,P29597 TRUE     
#>  4 O75791  Q13445        FALSE      O75791,Q13445 FALSE    
#>  5 O75791  Q6P4F2        FALSE      O75791,Q6P4F2 FALSE    
#>  6 O75791  Q8N9T8        FALSE      O75791,Q8N9T8 FALSE    
#>  7 O75791  Q9BYD3        FALSE      O75791,Q9BYD3 FALSE    
#>  8 O75791  Q9UMF0        FALSE      O75791,Q9UMF0 FALSE    
#>  9 P05362  O75791        TRUE       O75791,P05362 TRUE     
#> 10 P05362  P0C6A0        FALSE      P05362,P0C6A0 FALSE    
#> # … with 62 more rows
ppi_protein_graph <- protein_df_to_graph(shared_ppi_df)
ppi_protein_groups <- graph_to_protein_groups(ppi_protein_graph)
ppi_protein_groups %>% str()
#> tibble [5 × 2] (S3: tbl_df/tbl/data.frame)
#>  $ group_member:List of 5
#>   ..$ 1: chr [1:4] "O75791" "P05362" "P29597" "Q9UMF0"
#>   ..$ 5: chr [1:2] "Q8N9T8" "Q9BYD3"
#>   ..$ 2: chr "P0C6A0"
#>   ..$ 3: chr "Q13445"
#>   ..$ 4: chr "Q6P4F2"
#>   ..- attr(*, "dim")= int 5
#>   ..- attr(*, "dimnames")=List of 1
#>   .. ..$ : chr [1:5] "1" "5" "2" "3" ...
#>  $ group_size  : Named int [1:5] 4 2 1 1 1
#>   ..- attr(*, "names")= chr [1:5] "1" "5" "2" "3" ...
plot(ppi_protein_graph)

sessionInfo

sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] epigraphdb_0.2.2 igraph_1.2.6     glue_1.4.2       purrr_0.3.4     
#> [5] dplyr_1.0.6      magrittr_2.0.1  
#> 
#> loaded via a namespace (and not attached):
#>  [1] highr_0.9         pillar_1.6.1      compiler_4.1.0    tools_4.1.0      
#>  [5] digest_0.6.27     jsonlite_1.7.2    evaluate_0.14     memoise_2.0.0    
#>  [9] lifecycle_1.0.0   tibble_3.1.2      pkgconfig_2.0.3   rlang_0.4.11     
#> [13] cli_2.5.0         rstudioapi_0.13   curl_4.3.1        yaml_2.2.1       
#> [17] pkgdown_1.6.1     xfun_0.23         fastmap_1.1.0     stringr_1.4.0    
#> [21] httr_1.4.2        knitr_1.33        generics_0.1.0    desc_1.3.0       
#> [25] fs_1.5.0          vctrs_0.3.8       systemfonts_1.0.2 gtools_3.8.2     
#> [29] rprojroot_2.0.2   tidyselect_1.1.1  R6_2.5.0          textshaping_0.3.4
#> [33] fansi_0.5.0       rmarkdown_2.8     ps_1.6.0          htmltools_0.5.1.1
#> [37] ellipsis_0.3.2    ragg_1.1.2        utf8_1.2.1        stringi_1.6.2    
#> [41] cachem_1.0.5      crayon_1.4.1