diff --git a/DESCRIPTION b/DESCRIPTION index b1f0900..3adaa96 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,9 @@ Package: compareRhythms Type: Package Title: Approaches to identify Differential Rhythmicity -Version: 0.99.3 +Version: 1.0.0 Author: Bharath Ananthasubramaniam -Maintainer: Bharath Ananthasubramaniam +Maintainer: Bharath Ananthasubramaniam Description: This package identifies features that have altered rhythms between two conditions, i.e., are differentially rhythmic. Hypothesis testing and model selection pipelines classify differentially @@ -24,7 +24,8 @@ Imports: assertthat (>= 0.2.1), edgeR (>= 3.28.0), DESeq2 (>= 1.26.0), - SummarizedExperiment (>= 1.16.1) + SummarizedExperiment (>= 1.16.1), + lme4 (>= 1.1.29) Suggests: testthat (>= 2.1.0), statmod (>= 1.2.2), diff --git a/NAMESPACE b/NAMESPACE index cbfcbb9..e6b747f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,3 +1,4 @@ # Generated by roxygen2: do not edit by hand export(compareRhythms) +export(compareRhythms_cosinor) diff --git a/R/compareRhythms.R b/R/compareRhythms.R index f19eaa6..e3d08e0 100644 --- a/R/compareRhythms.R +++ b/R/compareRhythms.R @@ -14,7 +14,9 @@ #' selection, "dodr" for analysis using [DODR::dodr], "limma" for #' linear-modeling approach based on \pkg{limma}, "voom" for linear-modeling #' approach for RNA-Seq using [limma::voom], "deseq2" for RNA-seq analysis -#' using \pkg{DESeq2}, and "edger" for RNA-seq analysis using \pkg{edgeR}. +#' using \pkg{DESeq2}, "edger" for RNA-seq analysis using \pkg{edgeR}, and +#' "cosinor" for simple cosinor-based analysis both for independent samples +#' and repeated samples. #' @param period The period of rhythm being tested (default = 24) #' @param rhythm_fdr The false discovery cutoff for finding rhythmic time series #' (default = 0.05) @@ -35,6 +37,8 @@ #' different methods (default = TRUE). #' @param outliers Boolean specifying if weights must be computed for each #' sample to account for outliers. Only used by method = "voom". +#' @param longitudinal Boolean specifying if repeated samples from one +#' experimental unit. Only used by method = "cosinor". #' #' @return A *data.frame* with the names of the differentially rhythmic #' features, the category it is classified under and optionally the rhythm @@ -47,7 +51,8 @@ compareRhythms <- function(data, exp_design, lengths=NULL, method = "mod_sel", period=24, rhythm_fdr = 0.05, compare_fdr = 0.05, amp_cutoff = 0.5, criterion = "bic", schwarz_wt_cutoff = 0.6, - just_classify = TRUE, robust = TRUE, outliers = FALSE + just_classify = TRUE, robust = TRUE, outliers = FALSE, + longitudinal = FALSE ) { assertthat::assert_that( @@ -58,7 +63,7 @@ compareRhythms <- function(data, exp_design, lengths=NULL, assertthat::are_equal(ncol(data), nrow(exp_design)), assertthat::has_name(exp_design, c("time", "group")), assertthat::noNA(data), - method %in% c("mod_sel", "limma", "dodr", "voom", "deseq2", "edger"), + method %in% c("mod_sel", "limma", "dodr", "voom", "deseq2", "edger","cosinor"), is.factor(exp_design$group), is.numeric(exp_design$time), length(levels(exp_design$group)) == 2, @@ -76,16 +81,25 @@ compareRhythms <- function(data, exp_design, lengths=NULL, schwarz_wt_cutoff >=0 & schwarz_wt_cutoff <=1.0, assertthat::is.flag(just_classify), assertthat::is.flag(robust), - assertthat::is.flag(outliers) + assertthat::is.flag(outliers), + assertthat::is.flag(longitudinal) ) if (method %in% c("deseq", "edger") && !is.null(lengths)) { assertthat::assert_that(all(lengths>0), msg = "All transcript lengths are not positive") } + if (method == "cosinor" && longitudinal) { + assertthat::assert_that(assertthat::has_name(exp_design, "ID")) + if (assertthat::has_name(exp_design, "ID")) { + assertthat::assert_that(is.factor(exp_design$ID)) + } + } + if (assertthat::has_name(exp_design, "batch")) { assertthat::assert_that(is.factor(exp_design$batch)) } + switch (method, mod_sel = compareRhythms_model_select(data = data, exp_design = exp_design, @@ -126,6 +140,10 @@ compareRhythms <- function(data, exp_design, lengths=NULL, edger = compareRhythms_edgeR(counts = data, exp_design = exp_design, lengths = lengths, period = period, rhythm_fdr = rhythm_fdr, amp_cutoff = amp_cutoff, compare_fdr = compare_fdr, - just_classify = just_classify) + just_classify = just_classify), + cosinor = compareRhythms_cosinor(data = data, exp_design = exp_design, period = period, + rhythm_fdr = rhythm_fdr, compare_fdr = compare_fdr, + amp_cutoff = amp_cutoff, just_classify = just_classify, + longitudinal = longitudinal) ) } diff --git a/R/compareRhythms_cosinor.R b/R/compareRhythms_cosinor.R new file mode 100644 index 0000000..95c0bd8 --- /dev/null +++ b/R/compareRhythms_cosinor.R @@ -0,0 +1,133 @@ +#' Run differential rhythmicity analysis for microarray using limma +#' +#' @param eset A matrix of expression values with gene in the rows and samples in columns +#' @inheritParams compareRhythms +#' @keywords internal +#' @export + +compareRhythms_cosinor <- function(data, exp_design, period, rhythm_fdr, + compare_fdr, amp_cutoff, just_classify, longitudinal) { + + group_id <- base::levels(exp_design$group) + + exp_design <- base::cbind(exp_design, + inphase = cos(2 * pi * exp_design$time / period), + outphase = sin(2 * pi * exp_design$time / period)) + + lmer_control <- lme4::lmerControl(check.conv.singular = lme4::.makeCC(action = "ignore", tol = formals(lme4::isSingular)$tol)) + + if ("batch" %in% colnames(exp_design)) { + + if (longitudinal) { + fit <- lapply(1:nrow(data), + function(i) list(lme4::lmer(data[i,]~(1|ID) + group + group:inphase + group:outphase + batch, data = exp_design, REML=FALSE, control = lmer_control), + lme4::lmer(data[i,]~(1|ID) + group + inphase + outphase + batch, data = exp_design, REML=FALSE, control = lmer_control), + lme4::lmer(data[i,]~(1|ID) + group + batch, data = exp_design, REML=FALSE, control = lmer_control))) + } else { + fit <- lapply(1:nrow(data), + function(i) list(lm(data[i,]~0 + group + group:inphase + group:outphase + batch, data = exp_design), + lm(data[i,]~0 + group + inphase + outphase + batch, data = exp_design), + lm(data[i,]~0 + group + batch, data = exp_design))) + } + + fit_coeffs <- vapply(fit, function(f){ + coefficients <- if (longitudinal) lme4::fixef(f[[1]]) else coef(f[[1]]) + names(coefficients) <- gsub("group", "", names(coefficients)) + names(coefficients) <- gsub(":", "_", names(coefficients)) + return(coefficients) + }, FUN.VALUE = double(7L)) + } else { + + if (longitudinal) { + fit <- lapply(1:nrow(data), + function(i) list(lme4::lmer(data[i,]~(1|ID) + group + group:inphase + group:outphase, data = exp_design, REML=FALSE, control = lmer_control), + lme4::lmer(data[i,]~(1|ID) + group + inphase + outphase, data = exp_design, REML=FALSE, control = lmer_control), + lme4::lmer(data[i,]~(1|ID) + group, data = exp_design, REML=FALSE, control = lmer_control))) + } else { + fit <- lapply(1:nrow(data), + function(i) list(lm(data[i,]~0 + group + group:inphase + group:outphase, data = exp_design), + lm(data[i,]~0 + group + inphase + outphase, data = exp_design), + lm(data[i,]~0 + group, data = exp_design))) + } + + fit_coeffs <- vapply(fit, function(f){ + coefficients <- if (longitudinal) lme4::fixef(f[[1]]) else coef(f[[1]]) + names(coefficients) <- gsub("group", "", names(coefficients)) + names(coefficients) <- gsub(":", "_", names(coefficients)) + return(coefficients) + }, FUN.VALUE = double(6L)) + + } + + fit_coeffs <- t(fit_coeffs) + rownames(fit_coeffs) <- rownames(data) + + rhythmic_in_either <- vapply(fit, function(f) { + d <- anova(f[[3]], f[[1]], test=ifelse(longitudinal,"LRT", "F")) + ifelse(longitudinal, d$`Pr(>Chisq)`[2], d$`Pr(>F)`[2]) + }, FUN.VALUE = double(1L)) + + names(rhythmic_in_either) <- rownames(data) + + adj_pval <- p.adjust(rhythmic_in_either, method = "fdr") + + results <- compute_model_params(fit_coeffs, group_id, type = "coef") + + results <- data.frame(results) + + results$id <- rownames(results) + + rownames(results) <- NULL + + results$max_amp <- pmax(results[, paste0(group_id[1], "_amp")], + results[, paste0(group_id[2], "_amp")]) + + results$adj_p_val_A_or_B <- adj_pval + + results <- results[(results$adj_p_val_A_or_B < rhythm_fdr) & + (results$max_amp >= amp_cutoff), ] + + assertthat::assert_that(assertthat::not_empty(results), + msg = "Sorry no rhythmic genes in either dataset for the thresholds provided.") + + results$max_amp <- NULL + + diff_rhy_results <- vapply(fit, + function(f) { + d <- anova(f[[2]], f[[1]], test=ifelse(longitudinal,"LRT", "F")) + ifelse(longitudinal, d$`Pr(>Chisq)`[2], d$`Pr(>F)`[2]) + }, FUN.VALUE = double(1L)) + names(diff_rhy_results) <- rownames(data) + + diff_rhy_results <- diff_rhy_results[results$id] + + results$adj_p_val_DR <- stats::p.adjust(diff_rhy_results, + method = "BH") + results$diff_rhythmic <- results$adj_p_val_DR < compare_fdr + + results$rhythmic_in_A <- results[, paste0(group_id[1], "_amp")] > amp_cutoff + + results$rhythmic_in_B <- results[, paste0(group_id[2], "_amp")] > amp_cutoff + + results$category <- base::mapply(categorize, + results$rhythmic_in_A, + results$rhythmic_in_B, + results$diff_rhythmic) + + + main_cols <- c("id", "category", "rhythmic_in_A", "rhythmic_in_B", + "diff_rhythmic") + + results <- results[, c(main_cols, + base::setdiff(colnames(results), main_cols))] + + if (just_classify) { + results <- results[, main_cols] + } + + rownames(results) <- NULL + colnames(results) <- gsub("A", group_id[1], colnames(results)) + colnames(results) <- gsub("B", group_id[2], colnames(results)) + + return(results) +} diff --git a/man/compareRhythms.Rd b/man/compareRhythms.Rd index 4441920..7864f40 100644 --- a/man/compareRhythms.Rd +++ b/man/compareRhythms.Rd @@ -17,7 +17,8 @@ compareRhythms( schwarz_wt_cutoff = 0.6, just_classify = TRUE, robust = TRUE, - outliers = FALSE + outliers = FALSE, + longitudinal = FALSE ) } \arguments{ @@ -34,7 +35,9 @@ methods "deseq" and "edgeR".} selection, "dodr" for analysis using \link[DODR:dodr]{DODR::dodr}, "limma" for linear-modeling approach based on \pkg{limma}, "voom" for linear-modeling approach for RNA-Seq using \link[limma:voom]{limma::voom}, "deseq2" for RNA-seq analysis -using \pkg{DESeq2}, and "edger" for RNA-seq analysis using \pkg{edgeR}.} +using \pkg{DESeq2}, "edger" for RNA-seq analysis using \pkg{edgeR}, and +"cosinor" for simple cosinor-based analysis both for independent samples +and repeated samples.} \item{period}{The period of rhythm being tested (default = 24)} @@ -64,6 +67,9 @@ different methods (default = TRUE).} \item{outliers}{Boolean specifying if weights must be computed for each sample to account for outliers. Only used by method = "voom".} + +\item{longitudinal}{Boolean specifying if repeated samples from one +experimental unit. Only used by method = "cosinor".} } \value{ A \emph{data.frame} with the names of the differentially rhythmic