## ----setup, echo=FALSE, results="hide"----------------------------------------
knitr::opts_chunk$set(tidy=FALSE, cache=TRUE,
                      dev="png",
                      message=TRUE, error=FALSE, warning=TRUE)

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

## ----vignettes, eval=FALSE----------------------------------------------------
#  browseVignettes("DifferentialRegulation")

## ----citation-----------------------------------------------------------------
citation("DifferentialRegulation")

## ----load, message=FALSE------------------------------------------------------
library(DifferentialRegulation)

## ----specify_data-dir---------------------------------------------------------
data_dir = system.file("extdata", package = "DifferentialRegulation")

## ----specify_directories------------------------------------------------------
# specify samples ids:
sample_ids = paste0("organoid", c(1:3, 16:18))
# set directories of each sample input data (obtained via alevin-fry):
base_dir = file.path(data_dir, "alevin-fry", sample_ids)

# Note that alevin-fry needs to be run with `--use-mtx` option to store counts in a `quants_mat.mtx` file.
path_to_counts = file.path(base_dir,"/alevin/quants_mat.mtx")
path_to_cell_id = file.path(base_dir,"/alevin/quants_mat_rows.txt")
path_to_gene_id = file.path(base_dir,"/alevin/quants_mat_cols.txt")

## ----specify_directories_EC---------------------------------------------------
path_to_EC_counts = file.path(base_dir,"/alevin/geqc_counts.mtx")
path_to_EC = file.path(base_dir,"/alevin/gene_eqclass.txt.gz")

## ----load_USA_counts----------------------------------------------------------
sce = load_USA(path_to_counts,
               path_to_cell_id,
               path_to_gene_id,
               sample_ids)

sce

## ----cell-type----------------------------------------------------------------
path_to_DF = file.path(data_dir,"DF_cell_types.txt")
DF_cell_types = read.csv(path_to_DF, sep = "\t", header = TRUE)
matches = match(colnames(sce), DF_cell_types$cell_id)
sce$cell_type = DF_cell_types$cell_type[matches]
table(sce$cell_type)

## ----load_EC_counts-----------------------------------------------------------
EC_list = load_EC(path_to_EC_counts,
                  path_to_EC,
                  path_to_cell_id,
                  path_to_gene_id,
                  sample_ids)

## ----samples_design-----------------------------------------------------------
design = data.frame(sample = sample_ids,
                    group = c( rep("3 mon", 3), rep("6 mon", 3) ))
design

## ----compute_PB_counts--------------------------------------------------------
PB_counts = compute_PB_counts(sce = sce,
                              EC_list = EC_list,
                              design =  design,
                              sample_col_name = "sample",
                              group_col_name = "group",
                              sce_cluster_name = "cell_type")

## ----rm_EC_list---------------------------------------------------------------
rm(EC_list)

## ----sce-test-----------------------------------------------------------------
# sce-based test:
set.seed(169612)
results_USA = DifferentialRegulation(PB_counts,
                                     EC = FALSE)

## ----EC-test------------------------------------------------------------------
# EC-based test:
set.seed(169612)
results_EC = DifferentialRegulation(PB_counts)

## ----names-results------------------------------------------------------------
names(results_EC)

## ----visualize_gene_results---------------------------------------------------
head(results_EC$Differential_results, 3)

## ----visualize_US_results-----------------------------------------------------
head(results_EC$US_results, 3)

## ----visualize_USA_results----------------------------------------------------
head(results_EC$USA_results, 3)

## ----visualize_convergence_results--------------------------------------------
results_EC$Convergence_results

## ----plot_pi------------------------------------------------------------------
plot_pi(results_EC,
        type = "US",
        gene_id = results_EC$Differential_results$Gene_id[1],
        cluster_id = results_EC$Differential_results$Cluster_id[1])

## ----plot_pi_USA--------------------------------------------------------------
plot_pi(results_EC,
        type = "USA",
        gene_id = results_EC$Differential_results$Gene_id[1],
        cluster_id = results_EC$Differential_results$Cluster_id[1])

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

