## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    crop = NULL,
    # Related to https://stat.ethz.ch/pipermail/bioc-devel/
    # 2020-April/016656.html
    dpi = 96,
    fig.width = 6,
    fig.height = 4,
    fig.keep = "all",
    max.height = "300px"
)

local({
    hook_output <- knitr::knit_hooks$get("output")
    knitr::knit_hooks$set(output = function(x, options) {
        if (!is.null(options$max.height)) {
            options$attr.output <- c(
                options$attr.output,
                sprintf('style="max-height: %s;"', options$max.height)
            )
        }
        hook_output(x, options)
    })
})

## ----load-packages------------------------------------------------------------
library(TiDEomics)
library(org.Mm.eg.db)

## ----load-example-------------------------------------------------------------
data(tutorial_data, package = "TiDEomics")
data(tutorial_sample_info, package = "TiDEomics")

tutorial_data[1:5, 1:5]

## ----input-data---------------------------------------------------------------
tutorial_data_2 <- tutorial_data
exp1_values <- as.matrix(tutorial_data[, -1])
gene_sd <- apply(exp1_values, 1, sd)

set.seed(5678)
noise <- matrix(
    rnorm(length(exp1_values)),
    nrow = nrow(exp1_values)
)
tutorial_data_2[, -1] <- exp1_values + noise * 0.5 * gene_sd
tutorial_data_2$Feature[seq_len(50)] <- paste0(
    tutorial_data_2$Feature[seq_len(50)], "_Exp2"
)

## ----overlap-summary----------------------------------------------------------
shared_features <- intersect(tutorial_data$Feature, tutorial_data_2$Feature)
length(shared_features)

## ----combine-tide-------------------------------------------------------------
tide_results <- list()
tide_results[["Exp1"]] <- prepare_tide(
    data = tutorial_data, sample_ann = tutorial_sample_info,
    feature_threshold = 0, keep = "threshold", residual_threshold = 100)
tide_results[["Exp2"]] <- prepare_tide(
    data = tutorial_data_2, sample_ann = tutorial_sample_info,
    feature_threshold = 0, keep = "threshold", residual_threshold = 100)

combined <- combine_tide(tide_results, assay = 2)
names(combined)
combined$filter_summary

## ----trend, fig.width = 5-----------------------------------------------------
plot_trend(combined$se_list_filtered,
    features = combined$shared_features[1:4], assay = 2,
    multi_layout = "wrap", independent_y = FALSE, ncol = 2
)

plot_trend(combined$se_list_filtered,
    features = combined$shared_features[1:4], assay = 2,
    multi_layout = "overlay", independent_y = TRUE, ncol = 2
)

plot_trend(combined$se_list_filtered,
    features = combined$shared_features[1:4], assay = 2,
    multi_layout = "vertical", independent_y = FALSE
)

plot_trend(combined$se_list_filtered,
    features = combined$shared_features[1:4], assay = 2,
    multi_layout = "horizontal", independent_y = TRUE
)

## ----pca, fig.width = 8, fig.height = 3---------------------------------------
PC <- plot_pca(combined$se_list_filtered,
    assay = 2,
    shared_only = TRUE,
    plot_screeplot = TRUE,
    plot_loadings = FALSE,
    plot_morepc = FALSE
)
PC$plot

## ----umap, fig.width = 8, fig.height = 3--------------------------------------
umap_out <- plot_umap(combined$se_list_filtered,
    seed = 1234, assay = 2,
    shared_only = FALSE
)
umap_out$plot

## ----cormat, fig.width = 8, fig.height = 4------------------------------------
plot_cor_matrix(combined$se_list_filtered,
    method = "spearman",
    assay = 2,
    shared_only = FALSE,
    label_rep = TRUE, label_batch = TRUE,
    cellwidth = 1.5, cellheight = 1.5
)

## ----wgcna-prep, fig.width = 4------------------------------------------------
inp1 <- prepare_WGCNA(combined$se_list_filtered_shared$Exp1,
    assay = 2, powers = seq(1, 20),
    networkType = "signed", RsquaredCut = 0.8
)
inp1$plot

inp2 <- prepare_WGCNA(combined$se_list_filtered_shared$Exp2,
    assay = 2, powers = seq(1, 20),
    networkType = "signed", RsquaredCut = 0.8
)
inp2$plot

## ----consensus----------------------------------------------------------------
cons <- run_consensus_WGCNA(
    nets = list(Exp1 = inp1, Exp2 = inp2),
    powers = c(
        Exp1 = inp1$custom_powerEstimate,
        Exp2 = inp2$custom_powerEstimate
    ),
    minModuleSize = 10
)

## ----consensus-hubs-----------------------------------------------------------
hubs <- extract_consensus_hubs(cons, top_n = 1)
head(hubs)

## ----consensus-metrics--------------------------------------------------------
module_metrics <- summarise_module_metrics(cons)
module_metrics$Consensus

## ----module-table-------------------------------------------------------------
cons_module <- WGCNA_module(cons$network, exclude_grey = TRUE)

## ----modules-v, fig.width = 8, fig.height = 9---------------------------------
plot_modules_v(cons_module, combined$se_list_merged_shared,
    assay = 2, scale = TRUE,
    height_ratio = 2)

## ----modules-h, fig.width = 10, fig.height = 5--------------------------------
hub_colours <- c(
    "Consensus hub" = "#0072B2",
    "Exp1 hub" = "#D55E00",
    "Exp2 hub" = "#009E73"
)

plot_modules_h(cons_module, combined$se_list_merged_shared,
    assay = 2, scale = TRUE,
    mark_features = hubs,
    feature_colours = hub_colours,
    profile_scale = "shared", # or "dataset", "free"
    group_annotation = "both", # or "name", "colour", "none"
    legend_position = "right", # or "bottom"
    heatmap_width = 5, heatmap_height = 8,
    profile_width = 3
)

## ----modules-h-selected, fig.width = 10, fig.height = 5-----------------------
module_background <- WGCNA_module(cons$network, exclude_grey = FALSE)$Feature

go_modules <- enrichGO_list(cons_module,
    OrgDb = org.Mm.eg.db,
    universe = module_background,
    category = "BP",
    pvalueCutoff = 1, qvalueCutoff = 1) # keep all terms for visualisation

plot_modules_h(cons_module, combined$se_list_merged_shared,
    assay = 2, scale = TRUE,
    mark_features = hubs,
    feature_colours = hub_colours,
    enrich_list = go_modules$all,
    enrich_category = c("BP", "Features"),
    enrich_wrap_width = 4,
    profile_scale = "free", # or "dataset", "shared"
    group_annotation = "colour", # or "name", "both", "none"
    legend_position = "right",
    heatmap_width = 3, heatmap_height = 8,
    profile_width = 3,
    selected_modules = c(2, 5)
)

## ----upset-h, fig.width = 7, fig.height = 3.5---------------------------------
set_list <- lapply(combined$se_list_filtered, rownames)
background_features <- unique(unlist(set_list))

feature_sets <- list(
    Exp1_Exp2 = intersect(set_list$Exp1, set_list$Exp2),
    Exp1 = setdiff(set_list$Exp1, set_list$Exp2),
    Exp2 = setdiff(set_list$Exp2, set_list$Exp1)
)

go_all <- enrichGO_list(feature_sets,
    OrgDb = org.Mm.eg.db,
    universe = background_features,
    category = "BP",
    pvalueCutoff = 1, qvalueCutoff = 1) # keep all terms for visualisation

plot_upset_h(set_list,
    enrich_list = go_all$all,
    enrich_category = "BP",
    enrich_top_n = 2,
    enrich_wrap_width = 7
)

## ----upset-v, fig.width = 5, fig.height = 5-----------------------------------
bp_terms <- unique(as.data.frame(go_all$all$BP)$ID)
enrich_term_colours <- stats::setNames(
    rep(c("#0072B2", "#D55E00", "#009E73", "black"),
        length.out = length(bp_terms)),
    bp_terms
)

plot_upset_v(set_list,
    enrich_list = go_all$all,
    enrich_category = "BP",
    enrich_top_n = 2,
    term_wrap_width = 50,
    enrich_term_order = "intersection", # or "rank"
    enrich_term_colours = enrich_term_colours,
    include_all = FALSE # see clusterProfiler::dotplot `includeAll` usage
)

## ----scatter, fig.width = 5, fig.height = 4-----------------------------------
exp1_stat <- data.frame(
    Feature = rownames(combined$se_list_filtered$Exp1),
    Mean = rowMeans(
        SummarizedExperiment::assay(combined$se_list_filtered$Exp1, 2),
        na.rm = TRUE
    )
)
exp2_stat <- data.frame(
    Feature = rownames(combined$se_list_filtered$Exp2),
    Mean = rowMeans(
        SummarizedExperiment::assay(combined$se_list_filtered$Exp2, 2),
        na.rm = TRUE
    )
)

plot_scatter(exp1_stat, exp2_stat,
    col1 = "Mean", col2 = "Mean",
    name1 = "Exp1 mean abundance",
    name2 = "Exp2 mean abundance",
    colour_threshold = c(x = 1, y = 1),
    label_threshold = c(x = 2.5, y = 2.5)
)

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