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

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

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

## ----eval = FALSE-------------------------------------------------------------
# if (!requireNamespace("devtools", quietly = TRUE)) {
#     install.packages("devtools")
# }
# devtools::install_github("Fanding-Zhou/QRscore")

## -----------------------------------------------------------------------------
library(QRscore)
library(DESeq2)

## -----------------------------------------------------------------------------
data("example_dataset_raw_3000_genes")

## -----------------------------------------------------------------------------
bulk_sparse_mat <- example_dataset_raw_3000_genes$COUNTS
ages <- example_dataset_raw_3000_genes$METADATA$AGE
ages[ages %in% c("20-29", "30-39")] <- "20-39"
ages[ages %in% c("40-49", "50-59")] <- "40-59"
ages[ages %in% c("60-69", "70-79")] <- "60-79"

## ----message=FALSE, warning=FALSE---------------------------------------------
col_means <- colMeans(bulk_sparse_mat, na.rm = TRUE)
col_zeros <- colMeans(bulk_sparse_mat == 0, na.rm = TRUE)
# The following threshold can be modified
col_ids <- which(col_means > 5 & col_zeros < 0.2)

bulk_df <- bulk_sparse_mat[, col_ids]
expr_mat <- t(bulk_df)
coldata <- data.frame(age = ages)
dds <- DESeqDataSetFromMatrix(
    countData = expr_mat,
    colData = coldata,
    design = ~age
)
dds <- estimateSizeFactors(dds)
est_size_factors = dds$sizeFactor

## -----------------------------------------------------------------------------
print("Number of kept genes after prefiltering:")
print(length(col_ids))

## -----------------------------------------------------------------------------
## define age groups and subset dataset
kept_samples <- coldata$age %in% c("20-39", "60-79")
expr_mat_1 <- expr_mat[, kept_samples] # directly use raw counts
est_size_factors_1 <- est_size_factors [kept_samples]
coldata_1 <- coldata[kept_samples,]

## -----------------------------------------------------------------------------
results_2_sample <- QRscoreGenetest(expr_mat_1,
    coldata_1,
    size_factors = est_size_factors_1,
    pairwise_test = TRUE,
    pairwise_logFC = TRUE, test_mean = TRUE,
    test_dispersion = TRUE, num_cores = 2,
    approx = "asymptotic"
)

## -----------------------------------------------------------------------------
results_2_sample_mean <- results_2_sample$mean_test
results_2_sample_DEG <-
    results_2_sample_mean[
        results_2_sample_mean$QRscore_Mean_adj_p_value < 0.05,
    ]
head(results_2_sample_DEG)

## -----------------------------------------------------------------------------
results_2_sample_var <- results_2_sample$var_test
results_2_sample_DDG <-
    results_2_sample_var[results_2_sample_var$QRscore_Var_adj_p_value < 0.05, ]
head(results_2_sample_DDG)

## -----------------------------------------------------------------------------
results_3_sample <- QRscoreGenetest(
    expr_mat, coldata$age, size_factors = est_size_factors,
    pairwise_test = FALSE,
    pairwise_logFC = TRUE, test_mean = TRUE,
    test_dispersion = TRUE, num_cores = 2,
    approx = "asymptotic"
)

## -----------------------------------------------------------------------------
### DEGs
results_3_sample_mean <- results_3_sample$mean_test
results_3_sample_DEG <- results_3_sample_mean[
    results_3_sample_mean$QRscore_Mean_adj_p_value < 0.05,
]
head(results_3_sample_DEG)

## -----------------------------------------------------------------------------
### DDGs
results_3_sample_var <- results_3_sample$var_test
results_3_sample_DDG <-
    results_3_sample_var[results_3_sample_var$QRscore_Var_adj_p_value < 0.05, ]
head(results_3_sample_DDG)

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