## ----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)
    })
})

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

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

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

## ----create-input-------------------------------------------------------------
data_obj <- create_input(
    data = tutorial_data,
    sample_ann = tutorial_sample_info)

data_obj

## ----view-data----------------------------------------------------------------
assays(data_obj)[["orig"]] |> as.data.frame() |> utils::head()

## ----coldata------------------------------------------------------------------
colData(data_obj) # sample annotation

## ----palette------------------------------------------------------------------
custom_palette <- c(
    "untreated" = "#1b9e77", "IFNbeta" = "#d95f02",
    "IFNgamma" = "#7570b3", "LPS" = "#e7298a"
)
set_custom_palette(custom_palette)

## ----abundance----------------------------------------------------------------
plot_distribution(data_obj, facet_by = "Group")

## ----norm0--------------------------------------------------------------------
data_obj <- normalise_to_start(data_obj)

## ----split--------------------------------------------------------------------
data_obj_list <- split_groups(data_obj)

## ----merge-rep----------------------------------------------------------------
data_obj_merged_list <- merge_replicates(data_obj_list)

## ----merge-groups-------------------------------------------------------------
data_obj_merged <- merge_groups(data_obj_merged_list)

## ----cormat-------------------------------------------------------------------
plot_cor_matrix(data_obj,
    method = "spearman",
    label_rep = TRUE, label_batch = TRUE,
    cellwidth = 2, cellheight = 2
)

## ----pca----------------------------------------------------------------------
PC <- plot_pca(data_obj,
    # pc1 = 1, pc2 = 2, # default to plot PC1 and PC2
    plot_screeplot = TRUE,
    plot_loadings = FALSE,
    plot_morepc = TRUE,
    circle = FALSE,
    morepc = 1:5
)
PC$p_list

## ----pca-3D-------------------------------------------------------------------
plot_pca_3D(PC$pca, pcs = 1:3)

## ----pca-eigencor-------------------------------------------------------------
PCAtools::eigencorplot(PC$pca,
    metavars = c("Group", "Time"),
    components = paste0("PC", 1:5),
    col = colorRampPalette(c("#3C5488FF", "white", "#E64B35FF"))(100),
    colCorval = "black"
)

## ----umap---------------------------------------------------------------------
umap <- plot_umap(data_obj, seed = 1234)
umap$p_list

## ----pca-umap-group, fig.width = 8, fig.height = 2.5--------------------------
plot_pca_by_group(data_obj, circle = TRUE, arrow = TRUE, legend_pos = "top")
# also accepts a list: plot_pca_by_group(data_obj_list)

plot_umap_by_group(data_obj, seed = 1234, legend_pos = "top")
# also accepts a list: plot_umap_by_group(data_obj_list)

## ----de-between-time, fig.width = 8, fig.height = 2, warning = FALSE----------
DE_between_time_out <- DE_between_time(data_obj, assay = 1,
    filter = 1, trend = FALSE)

## ----de-between-time-all------------------------------------------------------
DE_between_time_out$all_list$IFNbeta$`t2-t0` |> utils::head()

## ----de-between-time-de-------------------------------------------------------
DE_between_time_out$de_list$IFNbeta$`t2-t0` |> utils::head()

## ----de-between-time-plot, fig.width = 7, fig.height = 2----------------------
# Filtered with thresholds in `DE_between_time()`
plot_DE_between_time(DE_between_time_out,
    fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3
)

# Re-filtering with new thresholds
plot_DE_between_time(DE_between_time_out,
    fontsize = 8, value = FALSE, nrow = 1, heatmap_width = 3,
    adjP_thres = 0.01, logFC_thres = 1
)

## ----de-between-group---------------------------------------------------------
DE_between_group_out <- DE_between_group(data_obj, assay = 2,
    filter = 1, trend = TRUE)

## ----de-between-group-plot, fig.width = 3, fig.height = 2, out.width = "70%"----
# Filtered with thresholds in `DE_between_group()`
plot_DE_between_group(DE_between_group_out)

# Re-filtering with new thresholds
plot_DE_between_group(DE_between_group_out, adjP_thres = 0.01, logFC_thres = 1)

## ----de-between-group-all-----------------------------------------------------
DE_between_group_out$all_list$`IFNgamma-untreated`$`24` |> utils::head()

## ----de-between-group-de------------------------------------------------------
DE_between_group_out$de_list$`IFNgamma-untreated` |> utils::head()

## ----volcano------------------------------------------------------------------
plot_volcano(DE_between_group_out,
    group1 = "untreated", group2 = "IFNgamma", time = 24,
    logFC_thres = 0.5, adjP_thres = 0.05, label = TRUE)

## ----property-----------------------------------------------------------------
data_obj_merged_list <- calc_feature_property(data_obj_merged_list,
    threshold = 0)
property_tb <- summarise_feature_property(data_obj_merged_list)

utils::head(property_tb)

## ----group-specific-----------------------------------------------------------
group_specific_features(property_tb, groups = c("untreated"),
    genename = FALSE, GO = FALSE
)

## ----impute-------------------------------------------------------------------
data_obj_merged_imp_list <- impute_groups(data_obj_merged_list)

## ----trendy-run---------------------------------------------------------------
set.seed(1234)
random_features <- sample(rownames(data_obj_merged_imp_list[[1]]), 50)

example_res_list <- run_Trendy(data_obj_merged_imp_list,
    feature = random_features,
    minExp = 0.5,
    maxK = 1,
    minNumInSeg = 2, meanCut = 0, NCores = 2
)

## ----trendy-plot, fig.height = 2.5, fig.width = 8-----------------------------
plot_segments(data_obj_merged_imp_list,
    example_res_list,
    feature = c("Slc25a51", "Aunip"), # example features
    ylab = "Log2(CPM + 1)"
)

## ----trendy-summary, fig.height = 2, fig.width = 4----------------------------
plot_breakpoints(example_res_list)

trendy_summary <- summarise_Trendy(example_res_list)
trendy_summary |> utils::head()

## ----trendy-list--------------------------------------------------------------
trendy_list <- extract_segment_trends(trendy_summary)
trendy_list$IFNbeta

## ----variance-decomposition, fig.height = 3-----------------------------------
# filter genes for variance decomposition:
# at least 50% values > 0 in at least 2 groups
decomp_filter_genes <- group_specific_features(property_tb,
    filter_ratio = 0.5,
    group_pct = 2 / 4,
    GO = FALSE, genename = FALSE
)$features

var_decomp <- decomp_variance(data_obj,
    features = decomp_filter_genes,
    assay = 1, core = 2
)

plot_variance(var_decomp, rank = "Time", top_n = 20)
plot_variance(var_decomp, rank = "Group", top_n = 20)

## ----wgcna-prepare------------------------------------------------------------
# Example filtering by residual variance < Q3
var_res_q3 <- stats::quantile(var_decomp$Residual, 0.75, na.rm = TRUE)
filter_wgcna <- var_decomp |> dplyr::filter(Residual < var_res_q3) |>
    dplyr::pull(Feature)

data_obj_wgcna <- data_obj[filter_wgcna, ]

## ----choose-power-------------------------------------------------------------
wgcna_input <- prepare_WGCNA(data_obj_wgcna, assay = 2,
    powers = seq(1, 20),
    networkType = "signed", RsquaredCut = 0.8
)
wgcna_input$plot
picked_power <- wgcna_input$powerEstimate
picked_power

## ----run-wgcna----------------------------------------------------------------
net <- run_WGCNA(wgcna_input,
    power = picked_power,
    # corType = "pearson", # other option is "bicor"
    numericLabels = TRUE
)

plot_WGCNA(net, fontsize = 8)

## ----extract-wgcna-modules----------------------------------------------------
gene_module <- WGCNA_module(net, exclude_grey = TRUE)

gene_module |>
    dplyr::group_by(Module) |>
    dplyr::summarise(n = dplyr::n())

## ----module-metrics-----------------------------------------------------------
summarise_module_metrics(net)

## ----plot-wgcna-modules, fig.width = 6, fig.height = 6------------------------
plot_modules_v(gene_module,
    data_obj_merged, scale = TRUE,
    ylabel = "Z-score of log2(CPM + 1)",
    height_ratio = 2
)

## ----go_rank, fig.width=5, fig.height=3---------------------------------------
gse_group <- enrichGO_rank(var_decomp,
    gene_rank_by = "Group",
    OrgDb = org.Mm.eg.db,
    keyType = "SYMBOL", category = "BP",
    go_rank_by = "p.adjust")

enrichplot::gseaplot2(gse_group[["BP"]], geneSetID = 1:3, base_size = 8)

## ----go-example, fig.width = 6------------------------------------------------
background_wgcna <- colnames(wgcna_input$data)

go_list <- enrichGO_list(
    gene_list = gene_module, OrgDb = org.Mm.eg.db,
    universe = background_wgcna,
    pvalueCutoff = 0.9, # get more results for demonstration
    qvalueCutoff = 0.9,
    simplify = FALSE,
    keyType = "SYMBOL"
)

plot_GO(go_list$all, plot_dotplot = TRUE,
    plot_emapplot = FALSE,
    plot_cnetplot = FALSE,
    showCategory_dotplot = 3)

## ----enrich-msigdb------------------------------------------------------------
# Mouse Hallmark gene sets
hallmark_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna,
    minGSSize = 5, category = "MH", species = "Mus musculus", db_species = "MM")

# Human C2 (chemical and genetic perturbations)
perturb_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna,
    category = "C2", subcategory = "CGP", species = "Mus musculus")

# Specific gene sets (category not required)
ifn_msigdb <- enrich_msigdb(gene_module, universe = background_wgcna,
    gene_sets = c("HALLMARK_INTERFERON_ALPHA_RESPONSE",
        "HALLMARK_INTERFERON_GAMMA_RESPONSE"),
    minGSSize = 5, species = "Mus musculus", db_species = "MM")

## ----plot-wgcna-modules-go, fig.width = 8, fig.height = 5---------------------
plot_modules_h(gene_module,
    data_obj_merged, scale = TRUE,
    ylabel = "Z-score of log2(CPM + 1)",
    enrich_list = go_list$all,
    enrich_category = "BP",
    heatmap_width = 6,
    heatmap_height = 8
)

## ----plot-wgcna-modules-multi, fig.width = 15, fig.height = 5-----------------
hub_features <- extract_hubs(net, top_n = 3)

plot_modules_h(gene_module,
    data_obj_merged, scale = TRUE,
    ylabel = "Z-score of log2(CPM + 1)",
    enrich_list = c(go_list$all, hallmark_msigdb),
    enrich_category = c("BP", "CC", "MH"),
    enrich_p_threshold = c(0.05, 0.05, 0.05),
    mark_features = hub_features,
    heatmap_width = 6,
    heatmap_height = 8
)

## ----plot-feature-------------------------------------------------------------
plot_trend(data_obj, assay = 1,
    features = c("Abtb1", "Dram1", "Ifi27", "Nufip1"),
    title = "Example features")

# Or pre-calculate mean and sd with calc_mean_sd()
table_mean_sd_orig <- calc_mean_sd(data_obj)$orig
plot_trend(table_mean_sd_orig,
    features = c("Abtb1", "Dram1", "Ifi27", "Nufip1"),
    title = "Example features")

## ----plot-residual------------------------------------------------------------
plot_trend(data_obj, assay = 1,
    features = var_decomp |>
        dplyr::arrange(Residual) |> utils::head(12) |> dplyr::pull(Feature),
    title = "Features with lowest residual variance"
)

plot_trend(data_obj, assay = 1,
    features = var_decomp |>
        dplyr::arrange(dplyr::desc(Residual)) |>
        utils::head(12) |> dplyr::pull(Feature),
    title = "Features with highest residual variance"
)

## ----prepare-tide-------------------------------------------------------------
tide <- prepare_tide(
    data = tutorial_data,
    sample_ann = tutorial_sample_info,
    filter_ratio = 0.5,
    min_groups = 2,
    keep = "below_quantile",
    residual_threshold = 0.75
)
tide$filter_summary

## ----tide-de------------------------------------------------------------------
tide$DE <- list(
    between_group = DE_between_group_out,
    between_time  = DE_between_time_out
)
tide$WGCNA <- net

tide$enrichment$GO_modules <- go_list
tide$enrichment$MSigDB <- hallmark_msigdb
tide$enrichment$GO_group <- gse_group

## ----isee, eval = FALSE-------------------------------------------------------
# iSEE::iSEE(data_obj)

## ----deedee-------------------------------------------------------------------
# Flatten nested results
de_flat <- flatten_DE(tide$DE)
enrich_flat <- flatten_enrich(tide$enrichment)

# Build DeeDeeExperiment
dde <- DeeDeeExperiment::DeeDeeExperiment(
    sce = tide$se,
    de_results = de_flat,
    enrich_results = enrich_flat
)

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

