## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE,
    fig.crop = FALSE)

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

## ----data---------------------------------------------------------------------
library(RhyME)
library(SummarizedExperiment)
library(BiocParallel)

example <- readRDS(system.file("extdata", "example_data.rds",
    package = "RhyME"))
se <- example$se
dim(se)
head(colData(se))

## ----fit----------------------------------------------------------------------
fit <- fitCosinor(se, subject_col = "subject", prior = "auto", BPPARAM = SerialParam())

## ----results------------------------------------------------------------------
result <- rowData(fit)
truth <- result$true_status == "Rhythmic"
selected_cpgs <- rownames(fit)[c(head(which(truth), 2),
    head(which(!truth), 2))]

columns <- c("true_status", "RhyME_pvalue", "RhyME_padj",
    "RhyME_amplitude", "RhyME_phase", "RhyME_mesor",
    "RhyME_converged")
result[selected_cpgs, columns]

## ----significant--------------------------------------------------------------
detected <- result$RhyME_padj < 0.05
significant <- result[detected, , drop = FALSE]
nrow(significant)

## ----summary------------------------------------------------------------------
metadata(fit)$RhyME$fitSummary
metadata(fit)$RhyME$ebPrior

## ----metacycle, results='hide'------------------------------------------------
Y <- data.frame(CycID = rownames(se), assay(se))
meta_fit <- MetaCycle::meta2d(
    infile = "simulated_beta", inDF = Y, filestyle = "csv",
    timepoints = se$time, minper = 24, maxper = 24,
    outputFile = FALSE)
meta_result <- meta_fit$meta

## ----align-metacycle, include=FALSE-------------------------------------------
meta_result <- meta_result[match(rownames(se), meta_result$CycID), ]

## ----method-comparison, echo=FALSE--------------------------------------------
pvalues <- data.frame(
    RhyME = result$RhyME_pvalue,
    JTK = meta_result$JTK_pvalue,
    LS = meta_result$LS_pvalue,
    MetaCycle = meta_result$meta2d_pvalue)
padj <- sapply(pvalues, p.adjust, method = "BH")
rownames(padj) <- rownames(se)
comparison <- do.call(rbind, lapply(colnames(padj), function(method) {
    called <- !is.na(padj[, method]) & padj[, method] < 0.05
    tp <- sum(called & truth)
    fp <- sum(called & !truth)
    data.frame(Method = method, TP = tp, FN = sum(truth) - tp,
        FP = fp, TN = sum(!truth) - fp,
        `Power (%)` = 100 * tp / sum(truth),
        `Observed FDP (%)` = if (sum(called)) 100 * fp / sum(called) else NA_real_,
        check.names = FALSE)
}))
knitr::kable(comparison, digits = c(0, 0, 0, 0, 0, 1, 2),
    row.names = FALSE, caption = "Known simulation truth versus calls at BH < 0.05.")

## ----plot, fig.width=8, fig.height=5.5----------------------------------------
plotCosinorFit(fit, cpgs = selected_cpgs, ncol = 2) +
    ggplot2::theme_gray()

## ----real-data----------------------------------------------------------------
beta_file <- system.file("extdata", "real_beta.csv.gz", package = "RhyME")
sample_file <- system.file("extdata", "real_samples.csv", package = "RhyME")
beta <- as.matrix(read.csv(beta_file, row.names = 1,
    check.names = FALSE))
sample_info <- read.csv(sample_file, row.names = 1)
sample_info <- sample_info[colnames(beta), , drop = FALSE]

real_se <- SummarizedExperiment(
    assays = list(beta = beta),
    colData = DataFrame(sample_info))
dim(real_se)
head(colData(real_se))

## ----minfi-input, eval=FALSE--------------------------------------------------
# se <- SummarizedExperiment(
#     assays = list(beta = minfi::getBeta(grSet)),
#     rowRanges = rowRanges(grSet),
#     colData = colData(grSet))

## ----fixed-default------------------------------------------------------------
se_independent <- example$se_independent
fit_fixed <- fitCosinor(se_independent, BPPARAM = SerialParam())
metadata(fit_fixed)$RhyME$ebPrior

## ----covariates---------------------------------------------------------------
subject_age <- setNames(seq(28, 56, by = 2), levels(se$subject))
se$age <- unname(subject_age[as.character(se$subject)])
i <- c(head(which(truth), 25), head(which(!truth), 25))
adjusted <- fitCosinor(se[i, ],
    time_col = "time", subject_col = "subject", covariates = "age",
    prior = "none", BPPARAM = SerialParam(), verbose = FALSE)

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