library(mvMAPIT)
 library(dplyr)
 library(ggplot2)
 library(ggrepel)
 library(tidyr)
 library(kableExtra)
 library(RcppAlgos)

mvmapit_hybrid <- mvmapit(
 t(data$genotype),
 t(data$trait),
 test = "hybrid"
)
fisher <- fishers_combined(mvmapit_hybrid$pvalues)

manhplot <- ggplot(fisher, aes(x = 1:nrow(fisher), y = -log10(p))) +
 geom_hline(yintercept = -log10(thresh), color = "grey40", linetype = "dashed") +
 geom_point(alpha = 0.75, color = "grey50") +
 geom_text_repel(aes(label=ifelse(p < thresh, as.character(id), '')))
 plot(manhplot)

thresh <- 0.05 / (nrow(fisher) / 2)
significant_snps <- fisher %>%
 filter(p < thresh) # Call only marginally significant SNPs
 
truth <- simulated_data$epistatic %>%
 ungroup() %>%
 mutate(discovered = (name %in% significant_snps$id)) %>%
 select(name, discovered) %>%
 unique()
 
significant_snps %>%
 mutate_if(is.numeric, ~(as.character(signif(., 3)))) %>%
 mutate(true_pos = id %in% truth$name) %>%
 kable(., linesep = "") %>%
 kable_material(c("striped"))

True Epistatic SNPs

To compare this list to the full list of causal epistatic SNPs of the simulations, refer to the following list. There are 5 causal SNPs. Of these 5 causal SNPs, 4 were succesfully discovered by mvMAPIT.

truth %>%
 kable(., linesep = "") %>%
 kable_material(c("striped"))

 # exhaustive search for p-values
pairs <- NULL
 if (nrow(significant_snps) > 1) {
 pairnames <- comboGeneral(significant_snps$id, 2)
 # Generate unique pairs of SNP names;
 # for length(names) = n, the result is a (n * (n-1)) x 2 matrix with one row corresponding to a pair
 for (k in seq_len(nrow(pairnames))) {
 fit <- lm(y ~ X[, pairnames[k, 1]]:X[, pairnames[k, 2]])
 p_value1 <- coefficients(summary(fit))[[1]][2, 4]
 p_value2 <- coefficients(summary(fit))[[2]][2, 4]
 tib <- dplyr::tibble(
 x = p_value1,
 y = p_value2,
 u = pairnames[k, 1],
 v = pairnames[k, 2]
 )
 pairs <- bind_rows(pairs, tib)
 }
}
 
 colnames(pairs) <- c(colnames(y), "var1", "var2")

Visualize Exhaustive Search Results

We plot the \(-\mathrm{log}_{10}(p)\) of the p-values for the regression coefficients as tile plot to highlight the identified interaction structure.

plotable <- pairs %>%
 pivot_longer(
 cols = starts_with("p_"),
 names_to = "trait",
 names_prefix = "trait ",
 values_to = "p",
 values_drop_na = TRUE
 ) %>%
 mutate(trait = case_when(
 trait == "p_01" ~ "Trait 1",
 trait == "p_02" ~ "Trait 2"))
tiles <- ggplot(data = plotable, aes(x=var1, y=var2, fill=-log10(p)))+
 geom_tile() +
 facet_wrap(~trait) +
 theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
 scale_fill_viridis_c()
 plot(tiles)

The only significant interactions after multiple testing correction are the interaction between snp00068 and snp_00665 as well as snp_00465 and snp_00917.

pairs %>%
 filter(p_01 < 0.05/15 | p_02 < 0.05/15) %>%
 kable(., linesep = "", digits = 14) %>%
 kable_material(c("striped"))

True epistataic structure

Compare the results of the exhaustive search to the true interaction structure. Notice that the only significant interactions in the exhaustive search are the two with the largest true effects.

true_interactions <- simulated_data$interactions %>%
 mutate(var1 = sprintf(group1, fmt = "snp_%05d")) %>%
 mutate(var2 = sprintf(group2, fmt = "snp_%05d")) %>%
 mutate(trait = case_when(
 trait == 1 ~ "Trait 1",
 trait == 2 ~ "Trait 2")) %>%
 select(-c("group1", "group2"))
X <- true_interactions[, c("var1", "var2")]
X <- t(apply(X, 1, sort))
true_interactions[,c("var1", "var2")] <- X
 
epistatic_pairnames <- comboGeneral(simulated_data$epistatic$name %>% unique(), 2)
true_pairs <- NULL
 for (k in seq_len(nrow(epistatic_pairnames))) {
 tib <- dplyr::tibble(var1 = epistatic_pairnames[k, 1],
 var2 = epistatic_pairnames[k, 2])
 true_pairs <- bind_rows(true_pairs, tib)
}
anti <- anti_join(true_pairs, true_interactions) %>%
 mutate(effect_size = 0)
 # Joining with `by = join_by(var1, var2)`
 
true_int_plot <- true_interactions %>%
 bind_rows(anti %>% mutate(trait = "Trait 1")) %>%
 bind_rows(anti %>% mutate(trait = "Trait 2"))
 
true_tiles <- ggplot(data = true_int_plot, aes(x=var1, y=var2, fill=effect_size)) +
 geom_tile() +
 facet_wrap(~trait) +
 theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
 scale_fill_gradient2(low = "blue", high = "red", mid = "white")
 plot(true_tiles)

While mvMAPIT does not identify the explicit partner, it still implicates more correct SNPs in this example. All true epistatic SNPs are listed in the following table.

true_interactions %>%
 kable(., linesep = "") %>%
 kable_material(c("striped"))