## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    echo = TRUE,
    tidy.opts = list(width.cutoff = 80),
    tidy = TRUE
)
# long links
adj_link <- paste0(
    "https://thomazbastiaanssen.github.io/anansi/",
    "articles/adjacency_matrices.html"
)
dif_link <- paste0(
    "https://thomazbastiaanssen.github.io/anansi/",
    "articles/differential_associations.html"
)

library(anansi)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("anansi")

## ----remotes-install, eval = FALSE--------------------------------------------
# install.packages("remotes")
# remotes::install_github("thomazbastiaanssen/anansi")

## ----load-libraries, eval=TRUE------------------------------------------------
# load anansi
library(anansi)

# load ggplot2 and ggforce to plot results
library(ggplot2)
library(ggforce)

# load example data + metadata from FMT Aging study
data(FMT_data)

## ----'ko-table-format', echo = FALSE------------------------------------------
cat("anansi expects samples to be rows and features to be columns:\n")
t(FMT_KOs[seq(4L),seq(10L)])

## ----'t-tables',  eval=TRUE---------------------------------------------------
# Transpose demo tables to input format
tableY <- t(FMT_metab)
tableX <- t(FMT_KOs)

## ----metadata-format, echo = FALSE--------------------------------------------
cat("anansi expects samples to be rows and features to be columns:\n")
`rownames<-`(FMT_metadata, NULL)[seq(10L),]

## ----weaveweb, eval=TRUE------------------------------------------------------
# Create AnansiWeb object from tables
web <- weaveWeb(
    x = "ko", y = "cpd", 
    link = kegg_link(),
    tableY = tableY,
    tableX = tableX,
    metadata = FMT_metadata
)

## ----bioc1--------------------------------------------------------------------
# Convert web to mae
mae <- asMAE(web)

## ----bioc2--------------------------------------------------------------------
web.mae <- weaveWeb(
    mae,
    tableY = "cpd", tableX = "ko",
    link = kegg_link()
)

out.mae <- anansi(
    web = web.mae,
    formula = ~Legend,
    adjust.method = "BH",
    verbose = TRUE
)

identical(web, web.mae)

## ----run-anansi, eval=TRUE----------------------------------------------------
# Run anansi on AnansiWeb object
anansi_out <- anansi(
    web = web,
    formula = ~Legend,
    adjust.method = "BH",
    verbose = TRUE
)

## ----show-output, eval=TRUE---------------------------------------------------
# General structure of anansi() output
head(anansi_out, c(5, 5))

# Let's take a look at the column names
colnames(anansi_out)

## ----bioc3--------------------------------------------------------------------
# Filter associations by q value
anansi_out <- anansi_out[anansi_out$full_q.values < 0.2, ]

# Visualize disjointed associations
plotAnansi(
    anansi_out,
    association.type = "disjointed",
    model.var = "Legend",
    signif.threshold = 0.05,
    fill_by = "group"
)

## ----r-web--------------------------------------------------------------------
rweb <- randomWeb(
    n_samples = 10,
    n_reps = 1,
    n_features_x = 8,
    n_features_y = 12,
    sparseness = 0.5
)

# The random object comes with some dummy metadata:
metadata(rweb)

# also see `?randomAnansi`

## ----pairwiseApply1-----------------------------------------------------------
# For each feature pair, was the value for x higher than the value for y?
pairwise_gt <- pairwiseApply(
    X = web,
    FUN = function(x, y) x > y, 
    USE.NAMES = TRUE
)

pairwise_gt[1:5, 1:5]

## ----pairwiseApply2-----------------------------------------------------------
# Run cor() on each pair of features
pairwise_cor <- pairwiseApply(
    X = web,
    FUN = function(x, y) cor(x, y), 
    USE.NAMES = FALSE
)

pairwise_cor

## ----plot-extra1--------------------------------------------------------------
# Use tidyr to wrangle the correlation r-values to a single column
library(tidyr)

anansiLong <- anansi_out |>
    pivot_longer(starts_with("All") | contains("FMT")) |>
    separate_wider_delim(
        name,
        delim = "_", names = c("cor_group", "param")
    ) |>
    pivot_wider(names_from = param, values_from = value)

# Now it's ready to be plugged into ggplot2, though let's clean up a bit more.
# Only consider interactions where the entire model fits well enough.
anansiLong <- anansiLong[anansiLong$full_q.values < 0.2, ]

ggplot(anansiLong) +

    # Define aesthetics
    aes(
        x = r.values,
        y = feature_X,
        fill = cor_group,
        alpha = disjointed_Legend_p.values < 0.05
    ) +

    # Make a vertical dashed red line at x = 0
    geom_vline(xintercept = 0, linetype = "dashed", colour = "red") +

    # Points show  raw correlation coefficients
    geom_point(shape = 21, size = 3) +

    # facet per compound
    ggforce::facet_col(~feature_Y, space = "free", scales = "free_y") +

    # fix the scales, labels, theme and other layout
    scale_y_discrete(limits = rev, position = "right") +
    scale_alpha_manual(
        values = c("TRUE" = 1, "FALSE" = 1 / 3),
        "Disjointed association\np < 0.05"
    ) +
    scale_fill_manual(
        values = c(
            "Young yFMT" = "#2166ac",
            "Aged oFMT" = "#b2182b",
            "Aged yFMT" = "#ef8a62",
            "All" = "gray"
        ),
        breaks = c("Young yFMT", "Aged oFMT", "Aged yFMT", "All"),
        name = "Treatment"
    ) +
    theme_bw() +
    ylab("") +
    xlab("Pearson's \u03c1")

## ----plot-extra2--------------------------------------------------------------
library(patchwork)

# Get a list containing data.frames for each feature pair.
feature_pairs <- getFeaturePairs(web)

# For each, prepare a ggplot object with those features on x & y axes.
plot_list <- lapply(
    feature_pairs,
    FUN = function(pair_i) {
        ggplot(pair_i) +
            aes(
                y = .data[[colnames(pair_i)[1L]]],
                x = .data[[colnames(pair_i)[2L]]]
            ) +
            # We'll use facet_grid for labels instead of labs().
            facet_grid(colnames(pair_i)[1L] ~ colnames(pair_i)[2L])
    }
)

# Feed the first six of our ggplots to patchwork, using wrap_plots()
wrap_plots(plot_list[seq_len(6L)], guides = "collect") &
    # Continue as with an ordinary ggplot, but use & instead of +
    geom_point(
        shape = 21, show.legend = FALSE,
        aes(fill = FMT_metadata$Legend)
    ) &
    geom_smooth(
        method = "lm", se = FALSE,
        aes(colour = FMT_metadata$Legend),
        show.legend = FALSE
    ) &

    # Appearance
    scale_fill_manual(
        aesthetics = c("colour", "fill"),
        values = c(
            "Young yFMT" = "#2166ac",
            "Aged oFMT" = "#b2182b",
            "Aged yFMT" = "#ef8a62",
            "All" = "gray"
        )
    ) &
    theme_bw() &
    labs(x = NULL, y = NULL)

## ----session-info-------------------------------------------------------------
sessionInfo()