## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  error    = FALSE,
  warning  = FALSE,
  eval     = TRUE,
  message  = FALSE
)

## ----create_dde, eval = FALSE-------------------------------------------------
# # do not run
# dde <- DeeDeeExperiment(
#   sce = sce, # sce is a SingleCellExperiment object
#   # de_results is a named list of de results
#   de_results = de_results,
#   # enrich_res is a named list of enrichment results
#   enrich_results = enrich_results
# )

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

## ----loadlib------------------------------------------------------------------
library("DeeDeeExperiment")

## ----loadlibraries------------------------------------------------------------
library("DeeDeeExperiment")

library("macrophage")
library("DESeq2")
library("edgeR")
library("DEFormats")

## ----load_data----------------------------------------------------------------
# load data
data(gse, package = "macrophage")

## ----dds_setup----------------------------------------------------------------
# set up design
dds_macrophage <- DESeqDataSet(gse, design = ~ line + condition)
rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15)
keep <- rowSums(counts(dds_macrophage) >= 10) >= 6

dds_macrophage <- dds_macrophage[keep, ]

## ----dds_subset---------------------------------------------------------------
# set seed for reproducibility
set.seed(42)
# sample randomly for 1k genes to speed up the processing
selected_genes <- sample(rownames(dds_macrophage), 1000)

dds_macrophage <- dds_macrophage[selected_genes, ]

## ----dds_macrophage-----------------------------------------------------------
dds_macrophage
# run DESeq
dds_macrophage <- DESeq(dds_macrophage)
# contrasts
resultsNames(dds_macrophage)

## ----extract_de_results-------------------------------------------------------
FDR <- 0.05
# IFNg_vs_naive
IFNg_vs_naive <- results(dds_macrophage,
                         name = "condition_IFNg_vs_naive",
                         lfcThreshold = 1,
                         alpha = FDR
)
IFNg_vs_naive <- lfcShrink(dds_macrophage,
                           coef = "condition_IFNg_vs_naive",
                           res = IFNg_vs_naive,
                           type = "apeglm"
)

# Salm_vs_naive
Salm_vs_naive <- results(dds_macrophage,
                         name = "condition_SL1344_vs_naive",
                         lfcThreshold = 1,
                         alpha = FDR
)
Salm_vs_naive <- lfcShrink(dds_macrophage,
                           coef = "condition_SL1344_vs_naive",
                           res = Salm_vs_naive,
                           type = "apeglm"
)

head(IFNg_vs_naive)
head(Salm_vs_naive)

## ----dds_to_dge---------------------------------------------------------------
# create DGE list
dge <- DEFormats::as.DGEList(dds_macrophage)

## ----limmarun-----------------------------------------------------------------
# normalize the counts
dge <- calcNormFactors(dge)

# create design for DE
design <- model.matrix(~ line + group, data = dge$samples)
# transform counts into logCPM
v <- voom(dge, design)

# fitting linear models using weighted least squares for each gene
fit <- lmFit(v, design)

## ----limma_contrasts----------------------------------------------------------
# available comparisons
colnames(design)
# setup comparisons
contrast_matrix <- makeContrasts(
  IFNg_vs_Naive = groupIFNg,
  Salm_vs_Naive = groupSL1344,
  levels = design
)

## ----ebayes-------------------------------------------------------------------
# apply contrast
fit2 <- contrasts.fit(fit, contrast_matrix)

# empirical Bayes moderation of standard errors
fit2 <- eBayes(fit2)
de_limma <- fit2 # MArrayLM object

## ----topTable-----------------------------------------------------------------
# show top 10 genes, first columns
topTable(de_limma, coef = "IFNg_vs_Naive", number = 10)[, 1:7] 

## ----glmfit-------------------------------------------------------------------
dge <- DEFormats::as.DGEList(dds_macrophage)

# estimate dispersion
dge <- estimateDisp(dge, design)
# perform likelihood ratio test
fit <- glmFit(dge, design)

## ----edgeR_contrasts----------------------------------------------------------
# available comparisons
colnames(design)
# setup comparisons
contrast_matrix <- makeContrasts(
  IFNg_vs_Naive = groupIFNg,
  Salm_vs_Naive = groupSL1344,
  levels = design
)

## ----dge_lrt------------------------------------------------------------------
# DGELRT objects
dge_lrt_IFNg_vs_naive <- glmLRT(fit,
                                contrast = contrast_matrix[, "IFNg_vs_Naive"])
dge_lrt_Salm_vs_naive <- glmLRT(fit,
                                contrast = contrast_matrix[, "Salm_vs_Naive"])

## ----topTag1------------------------------------------------------------------
# show top 10 genes, first few columns
topTags(dge_lrt_IFNg_vs_naive, n = 10)[, 1:7] 

## ----dge_exact----------------------------------------------------------------
# perform exact test
# exact test doesn't handle multi factor models, so we have to subset

# IFNg vs naive
keep_samples <- dge$samples$group %in% c("naive", "IFNg")
dge_sub <- dge[, keep_samples]
# droplevel
group <- droplevels(dge_sub$samples$group)
dge_sub$samples$group <- group

# recalculate normalization and dispersion
dge_sub <- calcNormFactors(dge_sub)
dge_sub <- estimateDisp(dge_sub, design = model.matrix(~group))
# exact test
# DGEExact object
dge_exact_IFNg_vs_naive <- exactTest(dge_sub, pair = c("naive", "IFNg")) 


# SL1344 vs naive
keep_samples <- dge$samples$group %in% c("naive", "SL1344")
dge_sub <- dge[, keep_samples]
# droplevel
group <- droplevels(dge_sub$samples$group)
dge_sub$samples$group <- group

# recalculate normalization and dispersion
dge_sub <- calcNormFactors(dge_sub)
dge_sub <- estimateDisp(dge_sub, design = model.matrix(~group))
# exact test
dge_exact_Salm_vs_naive <- exactTest(dge_sub, pair = c("naive", "SL1344"))

## ----topTag_2-----------------------------------------------------------------
topTags(dge_exact_Salm_vs_naive, n = 10)[, 1:7]

## ----load_FEAs----------------------------------------------------------------
# load Functional Enrichment Analysis results examples
data("topGO_results_list", package = "DeeDeeExperiment")
data("gost_res", package = "DeeDeeExperiment")
data("clusterPro_res", package = "DeeDeeExperiment")

## ----topGO--------------------------------------------------------------------
# ifng vs naive
head(topGO_results_list$ifng_vs_naive)

## ----gost---------------------------------------------------------------------
# salmonella vs naive
head(gost_res$result)

## ----clusterPro---------------------------------------------------------------
# ifng vs naive
head(clusterPro_res$ifng_vs_naive)

## ----fea_formats--------------------------------------------------------------
DeeDeeExperiment::supported_fea_formats()

## ----construct_dde------------------------------------------------------------
# initialize DeeDeeExperiment with dds object and DESeq results (as a named list)
dde <- DeeDeeExperiment(
  sce = dds_macrophage,
  de_results = list(
    IFNg_vs_naive = IFNg_vs_naive,
    Salm_vs_naive = Salm_vs_naive
  )
)

dde

## ----addDEA-------------------------------------------------------------------
# add a named list of results from edgeR
dde <- addDEA(dde, dea = list(
  dge_exact_IFNg_vs_naive = dge_exact_IFNg_vs_naive,
  dge_lrt_IFNg_vs_naive = dge_lrt_IFNg_vs_naive
))

# add results from limma as a MArrayLM object
dde <- addDEA(dde, dea = de_limma)

dde
# inspect the columns of the rowData
names(rowData(dde))

## ----addDEA_force-------------------------------------------------------------
dde <- addDEA(dde, dea = list(same_contrast = dge_lrt_Salm_vs_naive))

# overwrite results with the same name
# e.g. if the content of the same object has changed
dde <- addDEA(dde,
               dea = list(same_contrast = dge_exact_Salm_vs_naive),
               force = TRUE)

dde

## ----dde_with_limma_multi-----------------------------------------------------
dde_limma <- DeeDeeExperiment(sce = dds_macrophage,
                              de_results = de_limma)

dde_limma

## ----convert_limma_to_list----------------------------------------------------
dde_limma_list <- limma_list_for_dde(fit = de_limma,
                                     number = nrow(de_limma),
                                     sort.by = "none")

names(dde_limma_list)

## ----add_multi_limma----------------------------------------------------------
dde_limma <- addDEA(dde_limma,
                    dea = dde_limma_list)

dde_limma

## ----addFEA-------------------------------------------------------------------
# add FEA results as a named list
dde <- addFEA(dde,
               fea = list(IFNg_vs_naive = topGO_results_list$ifng_vs_naive))

# add FEA results as a single object
dde <- addFEA(dde, fea = gost_res$result)

# add FEA results and specify the FEA tool
dde <- addFEA(dde, fea = clusterPro_res, fea_tool = "clusterProfiler")

dde

## ----link_fea-----------------------------------------------------------------
dde <- linkDEAandFEA(dde,
                        dea_name = "Salm_vs_naive",
                        fea_name = c("salmonella_vs_naive", "gost_res$result")
)

## ----summary-linked-----------------------------------------------------------
summary(dde)

## ----addScenarioInfo----------------------------------------------------------
dde <- addScenarioInfo(dde,
                         dea_name = "IFNg_vs_naive",
                         info = "This results contains the output of a Differential Expression Analysis performed on data from the `macrophage` package, more precisely contrasting the counts from naive macrophage to those associated with IFNg."
)

## ----summary_mini-------------------------------------------------------------
# minimal summary
summary(dde)

## ----summary_fdr--------------------------------------------------------------
# specify FDR threshold for subsetting DE genes based on adjusted p-values
summary(dde, FDR = 0.01)

## ----summary_show_scenario----------------------------------------------------
# show contextual information, if available
summary(dde, show_scenario_info = TRUE)

## ----rename-------------------------------------------------------------------
# rename dea, one element
dde <- renameDEA(dde,
                  old_name = "de_limma",
                  new_name = "ifng_vs_naive_&_salm_vs_naive"
)

# multiple entries at once
dde <- renameDEA(dde,
                  old_name = c("dge_exact_IFNg_vs_naive", "dge_lrt_IFNg_vs_naive"),
                  new_name = c("exact_IFNg_vs_naive", "lrt_IFNg_vs_naive")
)

# rename fea
dde <- renameFEA(dde,
                  old_name = "gost_res$result",
                  new_name = "salm_vs_naive"
)

## ----remove-------------------------------------------------------------------
# removing dea
dde <- removeDEA(dde, c(
  "ifng_vs_naive_&_salm_vs_naive",
  "exact_IFNg_vs_naive",
  "dge_Salm_vs_naive"
))

# removing fea
dde <- removeFEA(dde, c("salmonella_vs_naive"))

dde

## ----names--------------------------------------------------------------------
# get DEA names
getDEANames(dde)

# get FEA names
getFEANames(dde)

## ----get_dea_res--------------------------------------------------------------
# access the 1st DEA if dea_name is not specified (default: minimal format)
getDEA(dde) |> head()

# access the 1st DEA, in original format
getDEA(dde, format = "original") |> head()

# access specific DEA by name, in original format
getDEA(dde,
    dea_name = "lrt_IFNg_vs_naive",
    format = "original"
) |> head()

# access specific DEA by name (default: minimal format)
getDEA(dde, dea_name = "lrt_IFNg_vs_naive") |> head()

## ----getDEAList---------------------------------------------------------------
# get dea results as a list, (default: minimal format)
lapply(getDEAList(dde), head)

## ----getDEAList-original, eval = FALSE----------------------------------------
# # get dea results as a list, in original format
# lapply(getDEAList(dde, format = "original"), head)

## ----getDEAInfo, eval = FALSE-------------------------------------------------
# # extra info
# getDEAInfo(dde)

## ----dea_info_specific--------------------------------------------------------
# retrieve specific information like the package name used for specific dea
dea_name <- "Salm_vs_naive"
getDEAInfo(dde)[[dea_name]][["package"]]

## ----get_fea_res--------------------------------------------------------------
# access the 1st FEA, (default: minimal format)
getFEA(dde) |> head()

# access the 1st FEA, in original format
getFEA(dde, format = "original") |> head()

# access specific FEA by name
getFEA(dde, fea_name = "ifng_vs_naive") |> head()

# access specific FEA by name, in original format
getFEA(dde, fea_name = "ifng_vs_naive", format = "original") |> head()

## ----getFEAList---------------------------------------------------------------
# get fea results as a list (default: minimal format)
lapply(getFEAList(dde), head)


# get all FEAs for this de_name, in original format
lapply(getFEAList(dde, dea_name = "IFNg_vs_naive", format = "original"), head)

## ----getFEAInfo, eval = FALSE-------------------------------------------------
# # extra info
# getFEAInfo(dde)

## ----fea_info_specific--------------------------------------------------------
# the tool used to perform the FEA
getFEAInfo(dde)[["ifng_vs_naive"]][["fe_tool"]]

## ----export_res, eval=FALSE---------------------------------------------------
# # export only DEAs
# out_dir <- "path/to/results_folder"
# 
# export_result_for_dde(dde,
#                       res_type = "dea",
#                       res_format = "minimal",
#                       output_dir = out_dir,
#                       file_name = "macrophage_dataset")
# 
# # export both DEAs and FEAs
# export_result_for_dde(dde,
#                       res_type = c("dea","fea"),
#                       res_format = "minimal",
#                       output_dir = out_dir,
#                       file_name = "macrophage_dataset")
# 
# # export DEAs, FEAs, and assays
# export_result_for_dde(dde,
#                       res_type = "all",
#                       res_format = "minimal",
#                       output_dir = out_dir,
#                       file_name = "macrophage_dataset",
#                       force = TRUE)
# 
# 

## ----save, eval = FALSE-------------------------------------------------------
# saveRDS(dde, "dde_macrophage.RDS")

## ----runiSEE, eval = FALSE----------------------------------------------------
# library("iSEE")
# iSEE::iSEE(dde)

## ----runGeneTonic, eval = FALSE-----------------------------------------------
# library("GeneTonic")
# my_dde <- dde
# 
# res_de <- getDEA(my_dde, dea_name = "IFNg_vs_naive", format = "original")
# res_enrich <- getFEA(my_dde, fea_name = "ifng_vs_naive")
# annotation_object <- mosdef::get_annotation_orgdb(dds_macrophage,
#                                                   "org.Hs.eg.db",
#                                                   "ENSEMBL")
# 
# GeneTonic::GeneTonic(
#   dds = dds_macrophage,
#   res_de = res_de,
#   res_enrich = res_enrich,
#   annotation_obj = annotation_object,
#   project_id = "GeneTonic_with_DeeDee"
# )

## ----sessioinfo---------------------------------------------------------------
sessionInfo()