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

## ----eval = FALSE-------------------------------------------------------------
#  if (!require("devtools", quietly = TRUE)) {
#      install.packages("devtools")
#  }
#  library(devtools)
#  install_github("wejlab/BatchQC")

## ----load---------------------------------------------------------------------
library(BatchQC)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  BatchQC()

## -----------------------------------------------------------------------------
data(signature_data)
data(batch_indicator)
se_object <- BatchQC::summarized_experiment(signature_data, batch_indicator)
SummarizedExperiment::assayNames(se_object) <- "log_intensity"
se_object$batch <- as.factor(se_object$batch)
se_object$condition <- as.factor(se_object$condition)

## -----------------------------------------------------------------------------
# CPM Normalization
se_object <- BatchQC::normalize_SE(se = se_object, method = "CPM",
    log_bool = FALSE, assay_to_normalize = "log_intensity",
    output_assay_name = "CPM_normalized_counts")

# CPM Normalization w/log
se_object <- BatchQC::normalize_SE(se = se_object, method = "CPM",
    log_bool = TRUE, assay_to_normalize = "log_intensity",
    output_assay_name = "CPM_log_normalized_counts")

# DESeq Normalization
se_object <- BatchQC::normalize_SE(se = se_object, method = "DESeq",
    log_bool = FALSE, assay_to_normalize = "log_intensity",
    output_assay_name = "DESeq_normalized_counts")

# DESeq Normalization w/log
se_object <- BatchQC::normalize_SE(se = se_object, method = "DESeq",
    log_bool = TRUE, assay_to_normalize = "log_intensity",
    output_assay_name = "DESeq_log_normalized_counts")

# log adjust
se_object <- BatchQC::normalize_SE(se = se_object, method = "none",
    log_bool = TRUE, assay_to_normalize = "log_intensity",
    output_assay_name = "log_normalized_counts")

## -----------------------------------------------------------------------------
# Combat correction
se_object <- BatchQC::batch_correct(se = se_object, method = "ComBat",
    assay_to_normalize = "log_intensity", batch = "batch",
    covar = c("condition"), output_assay_name = "Combat_corrected")

## -----------------------------------------------------------------------------
batch_design_tibble <- BatchQC::batch_design(se = se_object, batch = "batch", 
    covariate = "condition")
batch_design_tibble

## -----------------------------------------------------------------------------
confound_table <- BatchQC::confound_metrics(se = se_object, batch = "batch")
confound_table

## -----------------------------------------------------------------------------
pearson_cor_result <- BatchQC::std_pearson_corr_coef(batch_design_tibble)
pearson_cor_result

## -----------------------------------------------------------------------------
cramers_v_result <- BatchQC::cramers_v(batch_design_tibble)
cramers_v_result

## -----------------------------------------------------------------------------
ex_variation <- batchqc_explained_variation(se = se_object,
    batch = "batch", condition = "condition", assay_name = "log_intensity")
EV_boxplot <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[1]])
EV_boxplot
EV_boxplot_type_2 <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[2]])
EV_boxplot_type_2

## -----------------------------------------------------------------------------
ex_variation <- batchqc_explained_variation(se = se_object,
    batch = "batch", condition = "condition", assay_name = "log_intensity")
EV_table <- BatchQC::EV_table(batchqc_ev = ex_variation[[1]])
EV_table
EV_table_type_2 <- BatchQC::EV_table(batchqc_ev = ex_variation[[2]])
EV_table_type_2

## -----------------------------------------------------------------------------
ex_variation <- batchqc_explained_variation(se = se_object,
    batch = "batch", condition = "condition", assay_name = "Combat_corrected")
EV_boxplot <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[1]])
EV_boxplot
EV_boxplot_type_2 <- BatchQC::EV_plotter(batchqc_ev = ex_variation[[2]])
EV_boxplot_type_2
EV_table <- BatchQC::EV_table(batchqc_ev = ex_variation[[1]])
EV_table
EV_table_type_2 <- BatchQC::EV_table(batchqc_ev = ex_variation[[2]])
EV_table_type_2

## -----------------------------------------------------------------------------
heatmaps <- BatchQC::heatmap_plotter(se = se_object, assay = "log_intensity",
    nfeature = 38, annotation_column = c("batch", "condition"),
    log_option = "FALSE")
correlation_heatmap <- heatmaps$correlation_heatmap
correlation_heatmap

## -----------------------------------------------------------------------------
heatmaps <- BatchQC::heatmap_plotter(se = se_object, assay = "log_intensity",
    nfeature = 38, annotation_column = c("batch", "condition"),
    log_option = FALSE)
heatmap <- heatmaps$topn_heatmap
heatmap

## -----------------------------------------------------------------------------
dendrogram_plot <- BatchQC::dendrogram_plotter(se = se_object,
    assay = "log_intensity",
    batch_var = "batch",
    category_var = "condition")
dendrogram_plot$dendrogram

## -----------------------------------------------------------------------------
circular_dendrogram_plot <- BatchQC::dendrogram_plotter(
    se = se_object, assay = "log_intensity", batch_var = "batch",
    category_var = "condition")
circular_dendrogram_plot$circular_dendrogram

## -----------------------------------------------------------------------------
pca_plot <- BatchQC::PCA_plotter(se = se_object, nfeature = 20, color = "batch",
    shape = "condition", assays = c("log_intensity", "Combat_corrected"),
    xaxisPC = 1, yaxisPC = 2, log_option = FALSE)
pca_plot$plot
pca_plot$var_explained

## -----------------------------------------------------------------------------
differential_expression <- BatchQC::DE_analyze(se = se_object,
    method = "limma", batch = "batch", conditions = c("condition"),
    assay_to_analyze = "log_intensity")

## -----------------------------------------------------------------------------
names(differential_expression)

## -----------------------------------------------------------------------------
head(differential_expression[["batch2"]])

## -----------------------------------------------------------------------------
pval_plotter(differential_expression)
head(pval_summary(differential_expression))

## -----------------------------------------------------------------------------
value <- round((max(abs(
    differential_expression[[length(differential_expression)]][, 1]))
    + min(
        abs(differential_expression[[length(differential_expression)]][, 1])))
    / 2)
volcano_plot(DE_results = differential_expression[["batch2"]],
    pslider = 0.05,
    fcslider = value)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  file_location <- "location/to/save/object"
#  
#  saveRDS(object = se_object, file = file_location)

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()