1. New cosinor-based analysis for both independent and longitudinal

   data.
2. Incremented version.
3. Updated email address.
4. Updated vignette and master file call to cosinor analysis.

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 <bharath.[email protected]>
+Maintainer: Bharath Ananthasubramaniam <bharath.[email protected]>
 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),

NAMESPACE

@@ -1,3 +1,4 @@
 # Generated by roxygen2: do not edit by hand
 
 export(compareRhythms)
+export(compareRhythms_cosinor)

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)
   )
 }

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)
+}