---
title: "Detecting GxGxE interactions with case-parent triads using E-GADGETS"
author: 
- name: Michael Nodzenski
  affiliation: 
  - Department of Biostatistics, University of North Carolina, Chapel Hill, NC
  - Graduate Partnerships Program, National Institutes of Health, Bethesda, MD
  - Biostatistics and Computational Biology Branch, National Institute of Environmental Health Sciences, Research Triangle Park, NC
  email: michael.nodzenski@gmail.com
date: "December 28, 2022"
package: epistasisGA
output: 
  BiocStyle::html_document:
    toc_float: true
    fig_width: 5
bibliography: library.bib
vignette: >
  %\VignetteIndexEntry{E-GADGETS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}  
---

```{r setup, echo=FALSE, results="hide"}
knitr::opts_chunk$set(tidy = FALSE, cache = FALSE, dev = "png", message = FALSE,
                      error = FALSE, warning = FALSE)
```	

# Introduction 

This vignette describes how to implement the E-GADGETS method using the epistasisGA package. The E-GADGETS method is used to mine genetic data from case-parent triads for higher-order gene-by-environment interactions, in which the joint effect of a set of single nucleotide polymorphisms (SNPs) depends on candidate non-genetic factors. We here refer to non-genetic factors as 'environmental exposures', or just 'exposures', regardless of whether those factors are actually agents encountered in the environment. For example, disease severity (i.e., high, medium, low) could be considered an 'exposure'. Additionally, we here refer to gene-by-environment interactions that involve multiple SNPs as 'GxGxE' interactions, and those that also involve multiple exposures as 'GxGxExE' interactions. E-GADGETS is an extension of the previously described GADGETS method [@GADGET2020], and we advise users to consult the 'GADGETS' vignette prior to this one for a more detailed explanation of how GADGETS, and by extension, E-GADGETS, works. 

# Implementing E-GADGETS

## Load Data

We begin our demonstration of E-GADGETS by loading a simulated example of case-parent triad data. Note that E-GADGETS requires case-parent triads and, unlike GADGETS, does not accommodate disease-discordant siblings. In the example data, we load simulated maternal, paternal, and case genotypes, as well as the exposures. These data represent 24 SNPs from 1,000 families. Rows correspond to families, and columns represent SNP genotypes. Genotypes must be coded as 0, 1, or 2. The exposure is a binary factor with two levels (0, 1). In the input genotype data, SNPs 6, 12, and 18 are simulated to jointly interact with the exposure. 

```{r}
library(epistasisGA)
data("case.gxe")
case <- as.matrix(case.gxe)
data("dad.gxe")
dad <- as.matrix(dad.gxe)
data("mom.gxe")
mom <- as.matrix(mom.gxe)
data("exposure")
```

## Pre-process Data

The second step is to pre-process the data. Below, we default to the assumption that SNPs located on the same biological chromosome are in linkage, but users can more carefully tailor this argument based on individual circumstances if desired. For the example data, the SNPs are drawn from chromosomes 10-13, with the columns sorted by chromosome, and 6 SNPs per chromosome. We therefore construct a vector as follows indicating this assumed linkage structure:

```{r}
ld.block.vec <- rep(6, 4)
```

This vector indicates that the input genetic data has 4 distinct linkage blocks, with SNPs 1-6 in the first block, 7-12 in the second block, 13-18 in the third block, and 19-24 in the fourth block. Note the ordering of the columns in the input data must be consistent with the specified structure.

Now, we can execute pre-processing:

```{r}
pp.list <- preprocess.genetic.data(case, father.genetic.data = dad,
                                   mother.genetic.data = mom,
                                   ld.block.vec = ld.block.vec, 
                                   categorical.exposures = exposure)
```

Note that, above, the exposure data were input for the argument `categorical.exposures`. In doing so, the software will treat each input column as a factor variable and, ultimately, create dummy variables for each level. On the other hand, if the exposure of interest were continuous, users would need to specify the argument `continuous.exposures`, which will not dummy code the input. Strictly speaking, since the example exposure is binary, dummy coding would not be different from the existing coding, so we could have specified either argument.

Also note that E-GADGETS can, in principle, accept multiple exposures, but we have not fully tested that idea; specifically, if there are multiple exposures and at least one is continuous. If users nevertheless wish to input multiple exposures in which one or more exposures are continuous, they should be mindful that we have not tested the software in that context and interpret results with caution. With that said, if users are interested in simultaneous effects of both continuous and categorical exposures, the continuous data should be included for `continuous.exposures` and the categorical data should be specified for `categorical.exposures`. 

## Run E-GADGETS 

We now run E-GADGETS to nominate SNP-sets whose joint effects appear to depend on the exposure using the `run.gadgets` function. A more detailed discussion of the `run.gadgets` function and its arguments is available in the 'GADGETS' vignette. Like GADGETS, E-GADGETS requires the user to pre-specify the number of SNPs that may jointly interact with the exposure, controlled by the `chromosome.size` argument. We recommend running the algorithm for a range of sizes (2-5 or 2-6), but for this small example, we will only consider 3-4. 

```{r, message = FALSE}
set.seed(100)
run.gadgets(pp.list, n.chromosomes = 5, chromosome.size = 3, 
       results.dir = "size3_res", cluster.type = "interactive",
       registryargs = list(file.dir = "size3_reg", seed = 1300),
       n.islands = 8, island.cluster.size = 4, 
       n.migrations = 2)

run.gadgets(pp.list, n.chromosomes = 5, chromosome.size = 4, 
       results.dir = "size4_res", cluster.type = "interactive", 
       registryargs = list(file.dir = "size4_reg", seed = 1400),
       n.islands = 8, island.cluster.size = 4, 
       n.migrations = 2)
```


Next, we condense the sets of results using the function `combine.islands`. Note that in addition to the results directory path, the function requires as input a data.frame indicating the rsIDs (or a placeholder name), reference, and alternate alleles for each SNP in the data passed to `preprocess.genetic.data` as well as the list output by `preprocess.genetic.data`. 

```{r}
data(snp.annotations.mci)
size3.combined.res <- combine.islands("size3_res", snp.annotations.mci,
                                      pp.list)
size4.combined.res <- combine.islands("size4_res", snp.annotations.mci,
                                      pp.list)

```

The E-GADGETS fitness score involves a standardization step, where the standardization is based on the fitness scores from 10,000 random SNP-sets and permutations of the exposure initiated at algorithm outset. The relevant means and standard deviations are stored in the `results.dir` directory, in the `null.mean.sd.info.rds` file. That information will be required to carry out permutation-based inferential procedures, demonstrated later. 

We now examine the results:

```{r}
library(magrittr)
library(knitr)
library(kableExtra)
kable(head(size3.combined.res)) %>%
  kable_styling() %>%
  scroll_box(width = "750px")

```

We see that we have identified the correct SNPs with the simulated GxGxE effect. In this example, with just few input SNPs, E-GADGETS was able to identify only the correct SNP-set as risk-associated and no others. In real applications, we anticipate E-GADGETS will nominate multiple distinct SNP-sets for further follow-up study.  

The elements of the output are similar to those from GADGETS. The most useful components are simply the identities of the SNPs in the nominated set(s). Similar to GADGETS, E-GADGETS suggests a risk associated allele for each SNP in the set; that nomination process is described in the E-GADGETS paper. Specifically note that the 'diff.vec' columns report the fitted mean genotype difference vector for the exposure level that appears most strongly associated with allele transmissions, the details of which are available in the E-GADGETS paper. For this vignette, it will suffice to mention that a positive value for a 'diff.vec' column suggests the variant allele for the corresponding SNP is more strongly risk associated than the reference allele, while a negative value suggests the reverse. Unlike GADGETS, E-GADGETS does not try to determine whether risk associated alleles follow a recessive pattern of inheritance, and instead automatically assumes that carrying one or more copies of the nominated allele is risk associated. (That is why the "allele.copies" columns are always "1+" for E-GADGETS, while they may be "1+" or "2" for GADGETS.) The `fitness.score` column does not have a straightforward interpretation except that SNP-sets with higher fitness presumably are more likely to have GxGxE effects.

A difference from GADGETS output is that E-GADGETS also outputs columns that start with "parent_comp" and "transmission_comp". Those columns report the coefficients from the parental and transmission-based components of the fitness score, respectively. The two fitness score components are described in the E-GADGETS paper, and we refer readers to that manuscript for full details. In this example, our exposure has two levels, so the 'parent_comp.exp1.beta' reports the model coefficient from the parental component associated with the first exposure level (here the value of the first exposure level is '0') and the 'parent_comp.exp2.beta' reports the coefficient associated with the second exposure level (here the value of the second exposure level is '1'). Likewise, the 'transmission_comp.exp1.snp1.beta' column reports the coefficient from the transmission-based component associated with the first SNP in the SNP-set for the first level of the exposure, 'transmission_comp.exp1.snp2.beta' the second SNP in the SNP-set, etc. 


## Run Permutation-based Tests 

Like GADGETS, E-GADGETS has an associated global test of association. The test assesses the null hypothesis that, among top-scoring SNP-sets returned by E-GADGETS, none contain any GxE or GxGxE effects. The test is similar to that from GADGETS, except here, we shuffle the exposure, rather than randomizing case/complement-sibling labels, and then re-run E-GADGETS to generate a null distribution of fitness scores. Note that this test assumes the input SNPs are independent of the candidate exposure under the null. 
We begin by creating 4 data sets with the observed exposure randomly re-assigned:

```{r}
set.seed(1400)
permute.dataset(pp.list, 
                permutation.data.file.path = "perm_data",
                n.permutations = 4)

```

Now we re-run E-GADGETS on each permuted data set. Note that we pass the same standardization information to `run.gadgets` as we used for the observed data. 

```{r}
#pre-process permuted data 
preprocessed.lists <- lapply(seq_len(4), function(permutation){
  
  exp.perm.file <- file.path("perm_data", 
                             paste0("exposure.permute", permutation, ".rds"))
  exp.perm <- readRDS(exp.perm.file)
  preprocess.genetic.data(case, father.genetic.data = dad,
                          mother.genetic.data = mom,
                          ld.block.vec = ld.block.vec,
                          categorical.exposures = exp.perm)
})


#specify chromosome sizes
chrom.sizes <- 3:4

#specify a different seed for each permutation
seeds <- 4:7

#run GADGETS for each permutation and size 
perm.res <- lapply(chrom.sizes, function(chrom.size){
  
  # grab standardization info 
  orig.res.dir <- paste0("size", chrom.size, "_res")
  st.info <- readRDS(file.path(orig.res.dir, "null.mean.sd.info.rds"))
  st.means <- st.info$null.mean
  st.sd <- st.info$null.se
  
  lapply(seq_len(4), function(permutation){
  
    perm.data.list <- preprocessed.lists[[permutation]]
    seed.val <- chrom.size*seeds[permutation]
    res.dir <- paste0("perm", permutation, "_size", chrom.size, "_res")
    reg.dir <- paste0("perm", permutation, "_size", chrom.size, "_reg")
    run.gadgets(perm.data.list, n.chromosomes = 5, 
           chromosome.size = chrom.size,
           results.dir = res.dir, cluster.type = "interactive", 
           registryargs = list(file.dir = reg.dir, seed = seed.val),
           n.islands = 8, island.cluster.size = 4,
           n.migrations = 2, null.mean.vec = st.means, 
           null.sd.vec = st.sd)
    
  })
  
})

#condense the results 
perm.res.list <- lapply(chrom.sizes, function(chrom.size){
  
  lapply(seq_len(4), function(permutation){
  
    perm.data.list <- preprocessed.lists[[permutation]]
    res.dir <- paste0("perm", permutation, "_size", chrom.size, "_res")
    combine.islands(res.dir, snp.annotations.mci, perm.data.list)
    
  })
  
})

```

After the null distribution of fitness scores has been generated, the global test of association can be run with exactly the same commands as GADGETS. The only difference is that, for E-GADGETS, the null hypothesis being tested is that none of the input SNPs interact with the exposure (either via GxE or GxGxE interactions). Note here that we base the test on the top chromosome of each size (`n.top.scores = 1`), but we recommend the default (`n.top.scores = 10`) for real applications.

```{r}
# chromosome size 3 results

# function requires a list of vectors of
# permutation based fitness scores
chrom3.perm.fs <- lapply(perm.res.list[[1]], 
                         function(x) x$fitness.score)
chrom3.list <- list(observed.data = size3.combined.res$fitness.score,
                     permutation.list = chrom3.perm.fs)

# chromosome size 4 results
chrom4.perm.fs <- lapply(perm.res.list[[2]], 
                         function(x) x$fitness.score)
chrom4.list <- list(observed.data = size4.combined.res$fitness.score,
                     permutation.list = chrom4.perm.fs)

# list of results across chromosome sizes, with each list 
# element corresponding to a chromosome size
final.results <- list(chrom3.list, chrom4.list)

# run global test
global.test.res <- global.test(final.results, n.top.scores = 1)
global.test.res$pval

```

Aside from the global test, for a candidate SNP-set, we can also compute a p-value corresponding to the null hypothesis that none of the component SNPs interact with the exposure (via either GxE or GxGxE interactions). In this example, we choose the highest ranking E-GADGETS chromosome to test and thereafter run the test on the same data. Therefore, the 'p-value' generated by the software will not have uniform distribution under the null, and should not be used for hypothesis tests. We instead refer to it as an 'h-value'. If the test were instead run on data independent from that used by E-GADGETS, the p-value would allow for valid inference. Regardless, smaller h-values indicate less consistency with the null. 

We carry out the test as follows:

```{r}
top.snps <- as.vector(t(size3.combined.res[1, 1:3]))
d3.st.info <- readRDS(file.path("size3_res/null.mean.sd.info.rds"))
d3.st.means <- d3.st.info$null.mean
d3.st.sd <- d3.st.info$null.se
set.seed(10)
GxE.test.res <- GxE.test(top.snps, pp.list, 
                         null.mean.vec = d3.st.means, 
                         null.sd.vec = d3.st.sd)
GxE.test.res$pval
```

The test indicates the SNP-set's fitness score is inconsistent with what would be expected under the null. 

## Visualize Results

We visualize results from E-GADGETS using the same network-plotting approach as GADGETS, except larger SNP and SNP-pair scores now indicate evidence for involvement in GxE or GxGxE interactions for the implicated SNPs or pairs, rather than epistasis.

```{r, fig.width = 14, fig.height = 12}
# vector of 95th percentile of null fitness scores max
chrom.size.thresholds <- global.test.res$max.perm.95th.pctl

# chromosome size 3 threshold
d3.t <- chrom.size.thresholds[1]

# chromosome size 4 threshold
d4.t <- chrom.size.thresholds[2]

# create results list 
obs.res.list <- list(size3.combined.res[size3.combined.res$fitness.score >= 
                                          d3.t, ], 
                     size4.combined.res[size4.combined.res$fitness.score >= 
                                          d4.t, ])

# list of standardization information for each chromosome size
d4.st.info <- readRDS(file.path("size4_res/null.mean.sd.info.rds"))
d4.st.means <- d4.st.info$null.mean
d4.st.sd <- d4.st.info$null.se
null.means.list <- list(d3.st.means, d4.st.means)
null.sds.list <- list(d3.st.sd, d4.st.sd)

set.seed(10)
graphical.scores <- compute.graphical.scores(obs.res.list, pp.list, 
                                             null.mean.vec.list = null.means.list,
                                             null.sd.vec.list = null.sds.list)
network.plot(graphical.scores, pp.list, graph.area = 200,
             node.size = 40, vertex.label.cex = 2)

```


# Cleanup and sessionInfo() {.unnumbered}

```{r, results="hide"}
#remove all example directories 
chrom.sizes <- 3:4
perm.reg.dirs <- as.vector(outer(paste0("perm", 1:4), 
                                 paste0("_size", chrom.sizes, "_reg"), 
                                 paste0))
perm.res.dirs <- as.vector(outer(paste0("perm", 1:4), 
                                 paste0("_size", chrom.sizes, "_res"), 
                                 paste0))
lapply(c("size3_res", "size3_reg", "size4_res", "size4_reg",
         perm.reg.dirs, perm.res.dirs, 
         "perm_data"), unlink, recursive = TRUE)


```

```{r}
#session information 
sessionInfo()

```

# References {.unnumbered}
