| Title: | Unified Inference of Circadian Phase from Bulk and Pseudobulk Transcriptomes |
|---|---|
| Description: | Provides a single SummarizedExperiment-based interface to published circadian phase (time-of-sample-collection) inference methods for bulk and donor-level pseudobulk transcriptomes, along with a shared cross-validation benchmarking harness, cross-dataset transfer testing, and a calibrated uncertainty layer: bootstrap and Bayesian confidence intervals validated for real coverage, a fixed or dynamically-selected clock-gene quality flag, and consensus phase estimation across every bundled engine. Does not implement a new phase-inference algorithm; unifies existing published methods that are otherwise scattered across incompatible standalone implementations. |
| Authors: | Alperen Uysal [aut, cre] (ORCID: <https://orcid.org/0000-0003-3777-6411>) |
| Maintainer: | Alperen Uysal <[email protected]> |
| License: | Artistic-2.0 |
| Version: | 0.99.3 |
| Built: | 2026-10-08 12:16:23 UTC |
| Source: | https://github.com/BiocStaging/circaTime |
estimatePhase()
List the bundled reference models available to estimatePhase()
availableReferenceModels()availableReferenceModels()
This list is the exhaustive set of models estimatePhase() can use.
It is currently species = "Mmusculus" only, covering all 12 GSE54650
tissues: Molecular Timetable and
tauFisher, plus their consensus_accuracy/consensus_equal combinations, on
platform = "microarray", and Molecular Timetable only on
platform = "rnaseq". ZeitZeiger and TimeSignatR are wrapped engines but are
deliberately absent here: they are reachable through fitZeitZeiger()/
predictZeitZeiger(), fitTimeSignatR()/predictTimeSignatR(),
benchmarkPhase(), and transferPhase() rather than as pre-trained
estimatePhase() models.
A data.frame: species, tissue, method, platform, cv_mae
(that model's own held-out cross-validated circular MAE in hours; NA if
not measured), n_amp_genes, n_residual_draws. Rows whose method starts
with consensus_ are the combined models reachable via
estimatePhase(method = "consensus", consensus_weights = ...).
availableReferenceModels()availableReferenceModels()
Phase 3 Step 3.3. This does not replace bootstrapPhaseCI()/
circularQuantileInterval()'s empirical interval - it derives a genuinely
different, complementary object from the same residual pool: a von Mises
posterior predictive density for a sample's true phase, centred on the point
estimate, with concentration (circularConcentration()) measured directly
from that pool rather than assumed. circularQuantileInterval() answers "what
interval did the observed residuals fall in"; this answers "what is the
probability the true phase is any given value" - useful whenever a caller
needs more than an interval endpoint, e.g. P(|error| < 2h) for a specific
downstream decision threshold, or a smooth density to plot.
bayesianPhasePosterior( point, residualPool, level = 0.8, period = 24, nGrid = 121 )bayesianPhasePosterior( point, residualPool, level = 0.8, period = 24, nGrid = 121 )
point |
Single numeric point phase estimate (or a vector, one per
sample), in |
residualPool |
Numeric vector of real circular residuals (signed
deviations from truth, as stored in |
level |
Nominal credible level for the returned interval (default 0.8). |
period |
Length of the cycle (default 24). |
nGrid |
Number of points in the returned density grid per sample (default 121, i.e. every 12 minutes over a 24h period). |
A list with kappa (fitted von Mises concentration, shared across all
points), concentration (the underlying mean resultant length, [0, 1]),
credible_lo/credible_hi (numeric vectors, one per point), and density
(a list of data.frame(phase, density), one per point, each integrating to
1 over [0, period)).
Because the credible interval here is the equal-tailed interval of a
parametric (unimodal, symmetric) fit to the pool's concentration, while
circularQuantileInterval()'s interval is the empirical quantile interval
of the pool itself, the two intervals agree exactly only if the pool really is
von Mises-distributed. Real residual pools are approximately but not exactly
von Mises (see benchmarks/run_bayesian_ci_benchmark.R for a direct
comparison on real data) - the intended use is as a complementary posterior
summary, not a replacement calibration target for the already-validated
empirical interval.
pool <- circularDiff(rnorm(200, sd = 1.5), 0) # a tight, roughly von Mises pool post <- bayesianPhasePosterior(point = c(6, 18), residualPool = pool) post$kappa post$credible_lopool <- circularDiff(rnorm(200, sd = 1.5), 0) # a tight, roughly von Mises pool post <- bayesianPhasePosterior(point = c(6, 18), residualPool = pool) post$kappa post$credible_lo
Cross-validated benchmark of phase-inference engines against ground truth
benchmarkPhase( se, truth_col, methods = "timetable", assay_name = NULL, donor_col = NULL, nFolds = 5, period = 24, seed = NULL, engineArgs = list() )benchmarkPhase( se, truth_col, methods = "timetable", assay_name = NULL, donor_col = NULL, nFolds = 5, period = 24, seed = NULL, engineArgs = list() )
se |
A |
truth_col |
Name of the |
methods |
Character vector of engine names to benchmark: |
assay_name |
Name of the assay to use; defaults to the first assay. |
donor_col |
Name of the |
nFolds |
Number of cross-validation folds, split by donor (default 5). |
period |
Length of the cycle (default 24). Passed through to every engine
whose |
seed |
Optional integer seed for reproducible fold assignment; see
|
engineArgs |
Named list, keyed by method, of extra arguments to pass to
that engine's |
A list with two elements: summary, a data.frame with one row per
method (method, circular_mae, circular_cor, n, n_failed); and
residuals, a data.frame with one row per (method, held-out sample)
(method, sample, donor, fold, truth, pred, circular_error,
aligned, fold_failed).
If an engine's fit/predict errors on a given fold (e.g. because that
fold's training data happens to violate an engine-specific requirement, such
as tauFisher's rhythmicity-detection step needing every clock hour
represented), that fold's test samples get pred = NA, fold_failed = TRUE,
and a message() is emitted - the rest of that method's folds still run and
count normally (n/n_failed in summary report how many held-out samples
got a real prediction vs. how many were lost to a fold failure). This never
fabricates a result for a failed sample (AGENTS.md Section 0.3).
circular_mae is the primary accuracy number. For evenly-sampled time-course
ground truth (e.g. GSE54650), circular_cor can be numerically unstable and
even strongly negative for accurate predictions - see the "Known instability"
section of circularCor() before relying on it as a benchmark gate.
set.seed(1) nGenes <- 15; nSamples <- 24 time <- rep(seq(0, 23, length.out = 6), 4) theta <- time / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) mat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.3) })) rownames(mat) <- paste0("gene", seq_len(nGenes)) colnames(mat) <- paste0("s", seq_len(nSamples)) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) SummarizedExperiment::colData(se)$ZT <- time bm <- benchmarkPhase(se, truth_col = "ZT", methods = "timetable", nFolds = 4, seed = 1) bm$summaryset.seed(1) nGenes <- 15; nSamples <- 24 time <- rep(seq(0, 23, length.out = 6), 4) theta <- time / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) mat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.3) })) rownames(mat) <- paste0("gene", seq_len(nGenes)) colnames(mat) <- paste0("s", seq_len(nSamples)) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) SummarizedExperiment::colData(se)$ZT <- time bm <- benchmarkPhase(se, truth_col = "ZT", methods = "timetable", nFolds = 4, seed = 1) bm$summary
Refits the named engine on nBoot donor-resampled-with-replacement copies of
the training data and predicts testMat each time; the spread of those
bootstrap predictions gives each test sample a circular confidence interval
(via circularQuantileInterval()), while a single fit on the original,
non-resampled training data gives the reported point estimate. This wraps
whichever registered engine's fit/predict functions (the same
.circaTimeEngines registry benchmarkPhase() uses); it does not invent a new
phase-inference method (AGENTS.md 0.3).
bootstrapPhaseCI( method, trainMat, trainTime, testMat, donors = NULL, nBoot = 50, level = 0.8, seed = NULL, fitArgs = list(), predictArgs = list(), period = 24, residualPool = NULL )bootstrapPhaseCI( method, trainMat, trainTime, testMat, donors = NULL, nBoot = 50, level = 0.8, seed = NULL, fitArgs = list(), predictArgs = list(), period = 24, residualPool = NULL )
method |
Name of a registered engine. |
trainMat, trainTime, testMat
|
As in the engine's own |
donors |
Optional character/factor vector, one donor ID per training
column ( |
nBoot |
Number of bootstrap replicates (default 50). |
level |
Nominal coverage of the returned interval, in |
seed |
Optional integer seed; RNG state is saved/restored (AGENTS.md Section 10). |
fitArgs, predictArgs
|
Extra named arguments passed to the engine's
|
period |
Length of the cycle (default 24). |
residualPool |
Optional numeric vector of independent circular residuals
(each in |
A bootstrap replicate's fit/predict call is allowed to fail (e.g.
tauFisher's MetaCycle step erroring when a resample happens to drop every
replicate of one clock hour - the same failure mode benchmarkPhase() isolates
per-fold). Failed replicates are dropped from that sample's interval rather
than aborting the whole bootstrap (AGENTS.md 0.3: report the failure, do not
fabricate or crash).
A data.frame with one row per test sample: sample, phase (point
estimate from the original unresampled fit), phase_lo, phase_hi, level,
n_boot_ok, n_boot_failed. The raw per-replicate draws are attached as
attr(., "draws"), an ncol(testMat) x nBoot matrix (NA for failed
replicates), so calibrationCurve() can recompute intervals at other levels
without refitting anything.
residualPool mattersRefitting on resampled training data only captures how sensitive the fit is to
which training donors were used - it says nothing about the error a genuinely
new test sample's own biological noise or the model's irreducible bias
contributes, which on real data is usually the dominant source of error.
Verified directly on GSE54650 (mouse 12-tissue atlas): a training-resampling-only
bootstrap for the timetable engine gave a stated 80% interval that actually
covered truth only ~48% of the time - badly overconfident, exactly the failure
mode AGENTS.md Phase 3 warns about ("uncalibrated intervals are worse than
none"). Supplying residualPool - an independent sample of real circular
residuals (e.g. circularDiff(pred, truth) from a prior cross-validated
benchmark run on data the current testMat was not scored against) adds one
resampled residual to every bootstrap draw, folding that test-time noise back
into the interval. This is a standard prediction-interval technique (residual/
bagging bootstrap, as used e.g. for random forest prediction intervals), not a
new phase-inference method (AGENTS.md 0.3).
set.seed(1) nGenes <- 15; nSamples <- 20 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) ci <- bootstrapPhaseCI("timetable", trainMat, trainTime, trainMat[, 1:5, drop = FALSE], nBoot = 5, seed = 1 ) ciset.seed(1) nGenes <- 15; nSamples <- 20 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) ci <- bootstrapPhaseCI("timetable", trainMat, trainTime, trainMat[, 1:5, drop = FALSE], nBoot = 5, seed = 1 ) ci
Recomputes circularQuantileInterval() at each of levels from the raw
bootstrap draws already computed by bootstrapPhaseCI() (its "draws"
attribute), and reports what fraction of the time the interval actually
contains the known ground truth. This is the calibration check AGENTS.md
Phase 3 task 3 requires before a stated interval can be trusted at all:
"uncalibrated intervals are worse than none."
calibrationCurve(ci, truth, levels = c(0.5, 0.7, 0.8, 0.9, 0.95), period = 24)calibrationCurve(ci, truth, levels = c(0.5, 0.7, 0.8, 0.9, 0.95), period = 24)
ci |
A |
truth |
Numeric vector of ground-truth phases, same length and order as
|
levels |
Numeric vector of nominal coverage levels to evaluate (default
|
period |
Length of the cycle (default 24). |
A data.frame with one row per level: nominal, empirical
(observed coverage fraction), n (non-NA intervals evaluated).
set.seed(1) nGenes <- 15; nSamples <- 20 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) ci <- bootstrapPhaseCI("timetable", trainMat, trainTime, trainMat[, 1:5, drop = FALSE], nBoot = 5, seed = 1 ) calibrationCurve(ci, truth = trainTime[1:5], levels = c(0.5, 0.8))set.seed(1) nGenes <- 15; nSamples <- 20 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) ci <- bootstrapPhaseCI("timetable", trainMat, trainTime, trainMat[, 1:5, drop = FALSE], nBoot = 5, seed = 1 ) calibrationCurve(ci, truth = trainTime[1:5], levels = c(0.5, 0.8))
Unsupervised engines (e.g. CYCLOPS-family, CHIRAL) recover sample ordering around
the circle only up to an arbitrary rotation and reflection (AGENTS.md 7.2). This
finds the rotation offset and reflection (out of the two possible orientations)
that best aligns est to ref, then returns the aligned estimate.
circularAlign(est, ref, period = 24, allow_reflection = TRUE)circularAlign(est, ref, period = 24, allow_reflection = TRUE)
est |
Numeric vector of unsupervised phase estimates, in |
ref |
Numeric vector of reference/ground-truth phases, same length as |
period |
Length of the cycle (default 24). |
allow_reflection |
Logical; if |
The rotation for each orientation is estimated as the circular mean of the
pairwise differences to ref; this is the standard closed-form alignment used in
the CYCLOPS/CHIRAL literature and is a close practical approximation to the
rotation that minimises circular MAE.
A list with elements aligned (numeric vector, same length as est),
rotation (the applied additive offset), reflected (logical), and mae (the
circular MAE of aligned vs ref, over the non-NA fitting subset).
truth <- c(0, 4, 8, 12, 16, 20) est <- wrapPhase(-truth + 5) # reflected and rotated circularAlign(est, truth)$rotationtruth <- c(0, 4, 8, 12, 16, 20) est <- wrapPhase(-truth + 5) # reflected and rotated circularAlign(est, truth)$rotation
The standard circular-statistics measure of how tightly a sample of phases
clusters around its own circular mean: 1 for perfect agreement, 0 for phases
spread uniformly around the whole circle. Equivalent to 1 - circular variance. This is the primitive circularWeightedMean()'s "concentration"
attribute already computes internally for a weighted combination of engine
estimates; this function computes the same quantity directly for any raw
sample of phases (or phase deviations), which is what Phase 3 Step 3.3's
Bayesian posterior (see bayesianPhasePosterior()) needs applied to a
reference model's residual pool.
circularConcentration(x, period = 24, na.rm = TRUE)circularConcentration(x, period = 24, na.rm = TRUE)
x |
Numeric vector of phase values (or circular deviations, e.g. from
|
period |
Length of the cycle (default 24). |
na.rm |
Logical; drop |
A single numeric value in [0, 1]. NA_real_ if no usable input
remains.
circularConcentration(c(0, 0, 0)) # 1: no spread at all circularConcentration(c(0, 6, 12, 18)) # ~0: spread evenly around the circlecircularConcentration(c(0, 0, 0)) # 1: no spread at all circularConcentration(c(0, 6, 12, 18)) # ~0: spread evenly around the circle
The circular analogue of Pearson correlation. Must be used instead of cor()
whenever either input is a phase (AGENTS.md 7.1).
circularCor(x, y, period = 24)circularCor(x, y, period = 24)
x, y
|
Numeric vectors of phase values (same units, same |
period |
Length of the cycle (default 24). |
A single numeric value in [-1, 1].
This coefficient centres each variable on its own circular mean before
comparing them. When one variable's values are (near-)evenly spaced around the
entire circle - e.g. ground truth from a time course sampled every 2h across
a full 24h day, such as GSE54650 - that variable's resultant vector length is
~0, its circular mean becomes numerically arbitrary (driven by floating-point
noise, not signal), and the resulting correlation can come out strongly
negative even for near-perfect predictions. This is a property of the
Jammalamadaka-SenGupta coefficient itself, not a bug: verified directly on
GSE54650-shaped ground truth (24 evenly-spaced timepoints), a prediction set
with 0.2h mean circular error still yields a coefficient around -0.46.
circularMAE() does not share this failure mode and should be treated as the
primary accuracy metric for evenly-sampled benchmark data; treat circularCor
as a secondary/diagnostic number on such data, not a reliability gate.
circularCor(c(0, 6, 12, 18), c(1, 7, 13, 19))circularCor(c(0, 6, 12, 18), c(1, 7, 13, 19))
The parametric counterpart to circularQuantileInterval(): instead of taking
empirical quantiles of a finite set of bootstrap draws, this fits a von Mises
distribution - the circular analogue of the normal distribution, and the
standard choice for Bayesian circular models (Fisher 1993) - to the
concentration already measured from real data, and returns the analytic
credible interval of that distribution centred on mu. Because a von Mises
distribution is unimodal and symmetric, the equal-tailed interval computed
here is also the shortest (highest-density) interval at the same level - no
separate HDI computation is needed.
circularCredibleInterval(mu, kappa, level = 0.8, period = 24)circularCredibleInterval(mu, kappa, level = 0.8, period = 24)
mu |
Single numeric point phase estimate (or a vector, recycled per
sample), in |
kappa |
Von Mises concentration parameter (non-negative). See
|
level |
Nominal credible level (default 0.8, matching
|
period |
Length of the cycle (default 24). |
A list with lo, hi (interval bounds in [0, period), same length
as mu).
circularCredibleInterval(mu = 12, kappa = 5, level = 0.8)circularCredibleInterval(mu = 12, kappa = 5, level = 0.8)
Signed circular difference between two phases
circularDiff(pred, truth, period = 24)circularDiff(pred, truth, period = 24)
pred, truth
|
Numeric vectors of phase values (same units, same |
period |
Length of the cycle (default 24). |
Numeric vector of signed differences in (-period/2, period/2], i.e. the
shortest signed rotation from truth to pred.
circularDiff(0.5, 23.5) # 1, not -23circularDiff(0.5, 23.5) # 1, not -23
Vectorised version of the error term inside circularMAE(), useful for residual
analysis (e.g. benchmarkPhase() per-sample output).
circularError(pred, truth, period = 24)circularError(pred, truth, period = 24)
pred, truth
|
Numeric vectors of phase values in |
period |
Length of the cycle (default 24, for hours). |
Numeric vector of per-sample circular absolute errors.
circularError(c(0.5, 12), c(23.5, 0))circularError(c(0.5, 12), c(23.5, 0))
Tests whether truth lies on the arc walked from lo to hi in the
increasing direction (wrapping past period if needed) - i.e. the same
direction circularQuantileInterval() builds its interval in. A
zero-width interval (lo == hi) only covers an exact match, not the whole
circle.
circularIntervalCovers(truth, lo, hi, period = 24)circularIntervalCovers(truth, lo, hi, period = 24)
truth |
Numeric vector of ground-truth phase values. |
lo, hi
|
Numeric vectors of interval bounds, same length as |
period |
Length of the cycle (default 24). |
Logical vector, NA wherever truth, lo, or hi is NA.
circularIntervalCovers(23.9, 23, 1) # TRUE: 23.9 is on the 23 -> 1 arc circularIntervalCovers(12, 23, 1) # FALSEcircularIntervalCovers(23.9, 23, 1) # TRUE: 23.9 is on the 23 -> 1 arc circularIntervalCovers(12, 23, 1) # FALSE
The primary error metric for this package (AGENTS.md 7.1). Computed as
pmin(abs(d), period - abs(d)) where d = pred - truth, so that e.g. a
prediction of 23.5 against a truth of 0.5 counts as 1 hour of error, not 23.
circularMAE(pred, truth, period = 24, na.rm = TRUE)circularMAE(pred, truth, period = 24, na.rm = TRUE)
pred, truth
|
Numeric vectors of phase values in |
period |
Length of the cycle (default 24, for hours). |
na.rm |
Logical; if |
A single numeric value: the mean circular absolute error.
circularMAE(c(0.5, 12), c(23.5, 0))circularMAE(c(0.5, 12), c(23.5, 0))
Circular mean direction
circularMean(x, period = 24, na.rm = TRUE)circularMean(x, period = 24, na.rm = TRUE)
x |
Numeric vector of phase values. |
period |
Length of the cycle (default 24). |
na.rm |
Logical; drop |
A single numeric value in [0, period): the mean direction of x.
circularMean(c(23, 1)) # ~0, not 12circularMean(c(23, 1)) # ~0, not 12
Used by bootstrapPhaseCI() to turn a sample's bootstrap replicate
predictions into a phase_lo/phase_hi interval. Centres the draws on
their circular mean, takes the empirical (1-level)/2 and 1-(1-level)/2
quantiles of the signed circular deviation from that centre (via
circularDiff(), so a deviation is always in (-period/2, period/2] and
wraparound never inflates the spread), then maps those quantile deviations
back onto the circle. This assumes the draws are unimodal around their
circular mean, which holds for a bootstrap distribution of predictions
clustered near the true phase, but would not hold for, say, a bimodal
rotation/reflection-ambiguous set (AGENTS.md 7.2) - resolve that ambiguity
before calling this.
circularQuantileInterval(x, level = 0.8, period = 24)circularQuantileInterval(x, level = 0.8, period = 24)
x |
Numeric vector of phase draws, in |
level |
Nominal interval coverage, in |
period |
Length of the cycle (default 24). |
A list with center (circular mean of x), lo, hi (interval
bounds in [0, period)). All NA_real_ if fewer than 2 non-NA draws
remain.
circularQuantileInterval(c(23, 0, 1, 23.5, 0.5))circularQuantileInterval(c(23, 0, 1, 23.5, 0.5))
The weighted generalisation of circularMean(): each phase is converted to a
unit vector, the vectors are summed with the supplied weights, and the
direction of the resultant is returned. Used by estimatePhase()'s
method = "consensus" to combine several engines' phase estimates (AGENTS.md
Section 6 Step 3.1). circularWeightedMean(x, rep(1, length(x))) is identical
to circularMean(x).
circularWeightedMean(x, w, period = 24, na.rm = TRUE)circularWeightedMean(x, w, period = 24, na.rm = TRUE)
x |
Numeric vector of phase values. |
w |
Numeric vector of non-negative weights, same length as |
period |
Length of the cycle (default 24). |
na.rm |
Logical; drop pairs where |
The resultant vector's length is also returned (as attribute
"concentration", in [0, 1]): 1 means every input phase agreed exactly, 0
means they cancelled out completely. This is the standard circular measure of
agreement, and for a consensus estimate it is a directly interpretable
"did the engines agree?" signal rather than a hidden intermediate.
A single numeric value in [0, period), with attribute
"concentration" (numeric, [0, 1]). NA_real_ (concentration NA) if no
usable input remains or all weights are zero.
# two engines disagreeing by 2h, the second trusted twice as much circularWeightedMean(c(6, 8), c(1, 2)) # perfect agreement -> concentration 1 attr(circularWeightedMean(c(6, 6), c(1, 1)), "concentration")# two engines disagreeing by 2h, the second trusted twice as much circularWeightedMean(c(6, 8), c(1, 2)) # perfect agreement -> concentration 1 attr(circularWeightedMean(c(6, 6), c(1, 1)), "concentration")
Tissue-normalised amplitude of the AGENTS.md 7.7 core clock gene set, used to flag samples where the circadian clock is too damped for a phase estimate to be meaningful (AGENTS.md 7.4). Each matched gene is z-scored within its tissue group (so amplitude is comparable across tissues with different baseline expression/variance), then the z-scored block is reduced with PCA; amplitude is the Euclidean norm of a sample's first two PC scores. This mirrors the "eigengene" amplitude used in the CYCLOPS/circadian-transcriptomics literature (Anafi et al. 2017) for the same QC purpose: PC1/PC2 capture the near-circular manifold that a co-oscillating gene set traces out across a cohort, so their combined magnitude reflects how strongly a given sample's clock genes co-oscillate, independent of what phase they point to. This is a QC/amplitude summary, not a new phase-inference algorithm (AGENTS.md 0.3).
clockAmplitude(exprMat, tissue = NULL, minGenes = 4)clockAmplitude(exprMat, tissue = NULL, minGenes = 4)
exprMat |
Genes x samples numeric matrix, log-scale, with gene symbol
rownames (any case; matched case-insensitively against |
tissue |
Optional character/factor vector, one tissue label per sample
( |
minGenes |
Minimum number of matched core clock genes required (default 4); errors below this since amplitude from too few genes is not a meaningful summary. |
A numeric vector, ncol(exprMat) long, of non-negative amplitude
values, with attributes genesUsed (character vector of matched rownames),
varExplained (named numeric, fraction of variance captured by PC1+PC2 per
tissue group), and center/scale/rotation (named lists, one element per
tissue group: the scale() center/scale vectors and prcomp() rotation
matrix used for that group's PCA) - the latter three let a caller reproduce
this exact projection for a new sample against an already-fitted reference
group, which is what estimatePhase() does for small (often n=1) inference
calls where there aren't enough samples in the call itself to fit a fresh PCA.
NA for samples in a tissue group with fewer than 3 members (too few for
PCA/z-scoring to be meaningful).
AGENTS.md 7.4 expects this amplitude to be low "in aged tissue, tumours, and
critically ill donors," and Phase 3's gate requires showing flagged samples have
demonstrably worse phase error. On the two labelled real datasets tested so far
this effect was not statistically significant: GSE54650 (12-tissue healthy
mouse atlas, no aged/damped-clock samples; p=0.083) and GSE84511 (young vs aged
mouse epidermal stem cells, real ZT ground truth; p=0.20 for the derived
threshold, and mean amplitude barely differs by age group, p=0.74 - see
benchmarks/results/phase3_amplitude_threshold_GSE84511_timetable.csv and
notes/decisions.md). The likely reason, per GSE84511's source
paper (Solanas et al. 2017 Cell) and its companion liver study (Sato et
al. 2017 Cell): aging leaves core clock GENE expression population-average
"virtually unaltered" and instead reprograms which downstream genes oscillate -
so a core-clock-gene PCA amplitude may not be the right proxy for the kind of
age-related damping the literature actually reports, even though aged samples
in GSE84511 did have ~2.5x higher real cross-validated phase error than young
ones. Kept as documented, honest scope for now rather than silently dropped or
rounded up (AGENTS.md 0.3); revisit if a dataset or amplitude definition
surfaces that separates the two groups.
set.seed(1) genes <- c("Arntl", "Per1", "Per2", "Cry1", "Nr1d1", "Dbp") mat <- matrix(rnorm(6 * 20), nrow = 6, dimnames = list(genes, paste0("s", 1:20))) amp <- clockAmplitude(mat)set.seed(1) genes <- c("Arntl", "Per1", "Per2", "Cry1", "Nr1d1", "Dbp") mat <- matrix(rnorm(6 * 20), nrow = 6, dimnames = list(genes, paste0("s", 1:20))) amp <- clockAmplitude(mat)
The Phase 3 Step 3.2 replacement for the fixed 16-gene amplitude flag
(clockAmplitude()). Rather than scoring every dataset against the same
core-clock panel, this selects the informative gene panel from the input
dataset itself: each gene is fitted with a cosinor against the cohort's
(estimated) phases, and the panel is the set of genes whose fit beats a
phase-permuted null at a chosen FDR. A sample's coherence is then the fraction
of its own deviation-from-gene-mean that the cohort's fitted rhythm explains -
1 for a sample sitting exactly on the cohort's circadian manifold, near 0 (or
negative) for a sample whose clock genes are damped, desynchronised, or simply
inconsistent with the phase it was assigned.
clockCoherence( exprMat, phase, tissue = NULL, period = 24, nPerm = 200, fdr = 0.05, minPanel = 20, minSamples = 10, seed = 1 )clockCoherence( exprMat, phase, tissue = NULL, period = 24, nPerm = 200, fdr = 0.05, minPanel = 20, minSamples = 10, seed = 1 )
exprMat |
Genes x samples numeric matrix, log-scale, with gene rownames. |
phase |
Numeric vector of phases in hours, |
tissue |
Optional character/factor vector, one label per sample. If
supplied, genes are z-scored and the panel is selected within each tissue
group separately (as in |
period |
Rhythm period in hours (default 24). |
nPerm |
Number of phase permutations forming the null (default 200). |
fdr |
Benjamini-Hochberg FDR for panel selection (default 0.05). |
minPanel |
Minimum panel size (default 20). If fewer genes clear the null,
the top |
minSamples |
Minimum samples per group (default 10); smaller groups get
|
seed |
Seed for the permutation null, so the panel is reproducible. |
A numeric vector, ncol(exprMat) long, of coherence scores (1 = perfect
fit to the cohort rhythm, lower = worse; can be negative). Attributes:
panel (named list, selected gene symbols per group), panelSize (named
integer), relaxed (named logical, TRUE where the null-based selection
fell through to the top-minPanel fallback), r2Cutoff (named numeric, the
smallest R-squared admitted, NA when relaxed).
The legacy flag scores amplitude on 16 core clock genes. Its documented failure
(see clockAmplitude()) is that aging and similar insults leave core clock
gene expression close to unaltered while reprogramming which downstream genes
oscillate - so a fixed core-clock panel is measuring the wrong thing. Selecting
the panel per dataset lets those reprogrammed, dataset-specific oscillators
carry the QC signal.
phase is normally the estimated phase from estimatePhase(), not ground
truth, so a sample whose phase was estimated badly will fit the cohort cosinor
badly and score low. That is the intended behaviour (the flag is meant to catch
exactly those samples) but it does mean coherence is not an independent
measurement of data quality - it is a self-consistency score. Validation
against held-out true phase (benchmarks/run_dynamic_qc_benchmark.R) is
therefore the only claim made for it here; do not read it as a general
sample-quality metric. Supplying known sampling times as phase removes the
circularity and turns it into a pure rhythm-strength score.
Unlike the legacy flag, this needs a cohort: the panel and the fitted rhythm
are both estimated across samples. Below minSamples (default 10) it returns
all NA rather than a number computed from too little data, and
estimatePhase() falls back to the legacy flag in that case. This is a
deliberate design choice over storing a fixed dynamic basis in the bundled
reference models - a stored panel would be "dynamic" only with respect to the
training set, which is the very thing Step 3.2 sets out to remove.
clockAmplitude() for the legacy fixed-panel flag,
deriveAmplitudeThreshold() for turning either score into a cutoff.
set.seed(1) nGenes <- 200; nSamples <- 24 t <- seq(0, 23, length.out = nSamples) truePhase <- runif(nGenes, 0, 24) mat <- t(vapply(seq_len(nGenes), function(g) { 2 * cos(2 * pi * (t - truePhase[g]) / 24) + rnorm(nSamples, sd = 1) }, numeric(nSamples))) dimnames(mat) <- list(paste0("g", seq_len(nGenes)), paste0("s", seq_len(nSamples))) coh <- clockCoherence(mat, phase = t) summary(coh)set.seed(1) nGenes <- 200; nSamples <- 24 t <- seq(0, 23, length.out = nSamples) truePhase <- runif(nGenes, 0, 24) mat <- t(vapply(seq_len(nGenes), function(g) { 2 * cos(2 * pi * (t - truePhase[g]) / 24) + rnorm(nSamples, sd = 1) }, numeric(nSamples))) dimnames(mat) <- list(paste0("g", seq_len(nGenes)), paste0("s", seq_len(nSamples))) coh <- clockCoherence(mat, phase = t) summary(coh)
Returns the AGENTS.md 7.7 core clock gene set, one row per gene, with common
aliases and a source citation. Matching against expression data rownames is
always case-insensitive (see clockAmplitude()) - the mouse orthologs of every
gene in this set share the same letters as the human symbol, differing only in
case, so a single table serves both species without a separate mouse column.
coreClockGenes()coreClockGenes()
A data.frame with columns symbol (canonical human HGNC symbol),
aliases (semicolon-separated), source (citation).
head(coreClockGenes())head(coreClockGenes())
Finds the clockAmplitude() cutoff below which phase error is meaningfully
worse, from held-out labelled data - AGENTS.md 7.4/Phase 3 task 4 requires this
be derived empirically, not from intuition.
deriveAmplitudeThreshold(amp, error, minLowFrac = 0.05)deriveAmplitudeThreshold(amp, error, minLowFrac = 0.05)
amp |
Numeric vector of |
error |
Numeric vector of per-sample circular error (e.g. from
|
minLowFrac |
Minimum fraction of samples the |
Candidate thresholds are deciles of amp (10th-50th percentile). For each
candidate, samples split into low (amp <= candidate) and high
(amp > candidate) groups; the candidate is scored by a one-sided
Wilcoxon rank-sum test of error(low) > error(high), converted to a
common-language effect size (probability a random low-group error exceeds a
random high-group error, in [0, 1]), subject to the low group holding at
least minLowFrac of samples - too small a low group is not a reliable flag,
however well it happens to separate this particular dataset. The candidate
with the highest effect size among those satisfying the fraction constraint is
returned.
A list: threshold (chosen amplitude cutoff), effect_size
(common-language effect size at that cutoff), p_value (Wilcoxon test
p-value), low_mae/high_mae (mean error in each group), low_frac
(fraction of samples flagged low_amplitude).
set.seed(1) amp <- c(runif(15, 0, 1), runif(15, 2, 4)) error <- c(runif(15, 3, 5), runif(15, 0, 1)) deriveAmplitudeThreshold(amp, error)set.seed(1) amp <- c(runif(15, 0, 1), runif(15, 2, 4)) error <- c(runif(15, 3, 5), runif(15, 0, 1)) deriveAmplitudeThreshold(amp, error)
The package's main entry point (AGENTS.md Section 6). Looks up a pre-trained
PhaseReferenceModel for the requested (species, tissue, method, platform) -
see availableReferenceModels() for what is currently bundled - predicts phase
with that engine's own predict function, attaches a circular confidence
interval built from the reference model's real held-out residual pool (no live
bootstrap refit at call time - see PhaseReferenceModel), a core-clock
amplitude projected onto the reference cohort's PCA basis (works even for a
single-sample call, unlike a fresh clockAmplitude() fit), and a QC flag.
estimatePhase( se, method, tissue, species = "Mmusculus", platform = "microarray", consensus_weights = c("accuracy", "equal"), assay_name = NULL, level = 0.8, qc_method = c("legacy", "dynamic"), qc_threshold = 0, qc_min_samples = 10, map_aliases = TRUE, data_species = species, map_orthologs = !identical(data_species, species) )estimatePhase( se, method, tissue, species = "Mmusculus", platform = "microarray", consensus_weights = c("accuracy", "equal"), assay_name = NULL, level = 0.8, qc_method = c("legacy", "dynamic"), qc_threshold = 0, qc_min_samples = 10, map_aliases = TRUE, data_species = species, map_orthologs = !identical(data_species, species) )
se |
A |
method |
Name of a bundled engine for the requested tissue/species (see
|
tissue |
Required, no default (AGENTS.md 7.5) - a tissue code from
|
species |
Default |
platform |
Default |
consensus_weights |
Only used when |
assay_name |
Name of the assay to use. If |
level |
Nominal coverage of |
qc_method |
|
qc_threshold |
Only used when |
qc_min_samples |
Only used when |
map_aliases |
Logical; if |
data_species |
Species of |
map_orthologs |
Logical, default |
Per the Section 6 contract, this never hard-errors on a per-sample basis - a
sample that can't be scored gets phase = NA and an informative phase_flag.
Hard errors are reserved for malformed input: an unrecognised
species/tissue/method combination, or an se missing gene rownames.
se with colData columns added: phase, phase_lo, phase_hi
(empirical residual-pool interval), phase_credible_lo, phase_credible_hi
(the Phase 3 Step 3.3 von Mises credible interval built on the same
residual pool - see bayesianPhasePosterior(); metadata(se)$circaTime_bayesian
records the fitted concentration), clock_amp, phase_flag (factor:
"ok" / "low_amplitude" / "gene_coverage" / "failed"). clock_amp
holds whichever score qc_method produced (amplitude or coherence -
metadata(se)$circaTime_qc records which). For method = "consensus",
metadata(se)$circaTime_consensus additionally records the per-engine
estimates and the weights used (see the Consensus section).
Only species = "Mmusculus" is currently bundled, trained from GSE54650
(platform = "microarray") and GSE54651 (platform = "rnaseq", its RNA-seq
counterpart - see inst/scripts/train_reference_models.R) - calling with any
other species/tissue/method/platform combination is a hard error listing what
IS available, not a silent guess (AGENTS.md 7.5/0.3). platform is never
inferred from se itself; picking the wrong one does not error, it just scores
against a model trained on a different assay technology than your data. Gene matching is case-insensitive symbol matching,
the same convention clockAmplitude()/coreClockGenes() already use.
data_species/map_orthologs (AGENTS.md Section 4 Step 1.2) let se carry
genes in a different species' symbol convention than the bundled model - see
those arguments below - but this remains one-to-one symbol/ortholog matching,
not full Ensembl-ID mapping within a single species (see standardizeGeneIds()
for that). A larger multi-species/tissue reference cohort (the circaTimeData
ExperimentHub package, AGENTS.md Section 9 Phase 4 task 4) is future work, not
silently substituted here.
Alongside the empirical phase_lo/phase_hi interval (residual-pool
bootstrap quantiles), every call also
returns phase_credible_lo/phase_credible_hi: a von Mises credible interval
at the same level, fit from the same reference model's residual pool via
bayesianPhasePosterior() rather than a second, independent computation. The
two intervals answer related but distinct questions - the empirical one is
"what interval did the observed residuals actually span", the Bayesian one is
"what is a smooth probability distribution over the true phase" - and are
expected to be close but not identical (see ?bayesianPhasePosterior's
Relationship section for when and why they can diverge).
metadata(se)$circaTime_bayesian records the fitted kappa and
concentration; call bayesianPhasePosterior() directly for the full
posterior density rather than just the interval endpoints.
method = "consensus")Rather than making the user pick one engine, method = "consensus" runs every
engine bundled for that (species, tissue, platform) and combines their phase
estimates on the circle via circularWeightedMean() (AGENTS.md Section 6
Step 3.1). Two weighting schemes are available via consensus_weights:
"accuracy" (default) - tissue-aware: each engine is weighted by
1/cvMAE, where cvMAE is that engine's own held-out cross-validated
circular MAE on this specific tissue, measured during training and stored
on the model (see PhaseReferenceModel). An engine that is more
accurate on liver than on brainstem is therefore weighted differently in
each.
"equal" - a plain circular mean of the engines' estimates, with no
accuracy prior at all. Useful as the honest baseline for whether the
accuracy weighting is buying anything.
The weights are not hidden: metadata(se)$circaTime_consensus records the
engines used, the stored_weights (the trained prior) and the
effective_weights actually applied to this call, each engine's own
per-sample per_engine_phase, the per-sample concentration (how strongly
the engines agreed, in [0, 1]) and n_engines_used. An engine that fails
its gene-coverage check on your data has its weight set to zero and the
remainder renormalised, rather than contributing an unusable phase; the
consensus is flagged "gene_coverage" only if no engine had coverage.
What consensus is available over is limited by what is bundled: mouse microarray has both Molecular Timetable and tauFisher, so consensus there combines two engines. ZeitZeiger and TimeSignature are deliberately not bundled as pre-trained models (their fitted objects are two orders of magnitude larger for no measured accuracy benefit), and the mouse RNA-seq platform bundles Molecular Timetable only, so no consensus model exists for it - requesting one is an error listing what does exist, not a silent fallback to a single engine.
mat <- matrix(rnorm(200), nrow = 100, dimnames = list(paste0("g", 1:100), c("s1", "s2"))) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) # random genes won't overlap the bundled model's trained gene set - this shows # the "never hard-error on a sample" contract: phase_flag is set, not an error. se <- estimatePhase(se, method = "timetable", tissue = "Liv") SummarizedExperiment::colData(se)$phase_flagmat <- matrix(rnorm(200), nrow = 100, dimnames = list(paste0("g", 1:100), c("s1", "s2"))) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) # random genes won't overlap the bundled model's trained gene set - this shows # the "never hard-error on a sample" contract: phase_flag is set, not an error. se <- estimatePhase(se, method = "timetable", tissue = "Liv") SummarizedExperiment::colData(se)$phase_flag
Fit a Molecular Timetable model from labelled training data
fitMolecularTimetable( trainMat, trainTime, period = 24, nGenes = 100, minCor = 0.5 )fitMolecularTimetable( trainMat, trainTime, period = 24, nGenes = 100, minCor = 0.5 )
trainMat |
Numeric matrix of expression, genes x samples (log-scale expression is assumed, as for the harmonic-regression fit to be meaningful). |
trainTime |
Numeric vector of known circadian phase (same units as
|
period |
Length of the cycle (default 24, for hours). |
nGenes |
Maximum number of time-indicating genes to keep, ranked by fit quality (Pearson correlation between fitted and observed expression). Default 100, following the general magnitude used in the original method. |
minCor |
Minimum per-gene fit correlation required to be eligible as a
time-indicating gene, evaluated before the |
An object of class "molecularTimetableFit": a list with the selected
gene names (genes), their harmonic-fit coefficients (mesor, b, c), fit
correlation (r), and period.
set.seed(1) nGenes <- 20; nSamples <- 24 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(nSamples, sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitMolecularTimetable(trainMat, trainTime, nGenes = 10) length(fit$genes)set.seed(1) nGenes <- 20; nSamples <- 24 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(nSamples, sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitMolecularTimetable(trainMat, trainTime, nGenes = 10) length(fit$genes)
Calls the authors' own pipeline (tauFisher::run_metacycle() +
find_periodic_genes() + train_tauFisher()) rather than reimplementing it.
run_metacycle() writes intermediate files; these are written to and cleaned up
from a temporary directory, never the working directory.
fitTauFisher( trainMat, trainTime, method = c("JTK", "LS"), thres = 1.01, nrep = 1, numbasis = 5, seed = NULL )fitTauFisher( trainMat, trainTime, method = c("JTK", "LS"), thres = 1.01, nrep = 1, numbasis = 5, seed = NULL )
trainMat |
Numeric matrix of expression, genes x samples (linear or
log-scale; |
trainTime |
Numeric vector of known circadian phase in |
method |
Character vector of MetaCycle method(s) to detect rhythmic genes,
passed to |
thres |
p-value threshold passed to |
nrep |
Number of training replicates, passed to |
numbasis |
Number of Fourier basis functions for the FDA smoothing step,
passed to |
seed |
Optional integer seed. |
An object of class "tauFisherFit": a list with the chosen periodic genes
(Genes), the PCA embedding (PCA_embedding) and trained multinomial model
(Model) returned by train_tauFisher(), which MetaCycle method supplied the
genes (method), and period (always 24).
if (requireNamespace("tauFisher", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 100 + 40 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitTauFisher(trainMat, time, seed = 1) length(fit$Genes) }if (requireNamespace("tauFisher", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 100 + 40 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitTauFisher(trainMat, time, seed = 1) length(fit$Genes) }
Calls the authors' own TimeSignatR::trainTimeStamp() rather than reimplementing
it. Unlike fitMolecularTimetable()/fitTauFisher(), TimeSignatR does no gene
selection of its own - the fitted cv.glmnet object's coefficient matrix spans
every gene in trainMat, with the elastic-net penalty (not a pre-filtering step)
doing the shrinkage - so predictTimeSignatR() requires every one of those genes
to be present in testMat, same all-or-nothing rule as predictTauFisher().
fitTimeSignatR( trainMat, trainTime, subjIDs = NULL, period = 24, alpha = 0.5, recalib = FALSE, seed = NULL, ... )fitTimeSignatR( trainMat, trainTime, subjIDs = NULL, period = 24, alpha = 0.5, recalib = FALSE, seed = NULL, ... )
trainMat |
Numeric matrix of expression, genes x samples. |
trainTime |
Numeric vector of known circadian phase in |
subjIDs |
Character/factor vector of subject/donor IDs, one per column of
|
period |
Length of the cycle (default 24, for hours). |
alpha |
Elastic-net mixing parameter, passed to |
recalib |
Logical, passed to |
seed |
Optional integer seed. |
... |
Additional arguments passed to |
An object of class "timeSignatRFit": a list with the training gene
names (genes), the raw trainTimeStamp() result (timestamp, containing the
fitted cv.glmnet object predTimeStamp() needs), and period.
if (requireNamespace("TimeSignatR", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) colnames(trainMat) <- paste0("s", seq_len(length(time))) fit <- fitTimeSignatR(trainMat, time, seed = 1) length(fit$genes) }if (requireNamespace("TimeSignatR", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) colnames(trainMat) <- paste0("s", seq_len(length(time))) fit <- fitTimeSignatR(trainMat, time, seed = 1) length(fit$genes) }
Calls the authors' own pipeline (zeitzeiger::zeitzeigerFit() +
zeitzeiger::zeitzeigerSpc()) rather than reimplementing it. Unlike
fitTauFisher(), zeitzeiger has no separate "trained model" object - the fitted
periodic splines are re-derived from xTrain/timeTrain at prediction time
(zeitzeiger::zeitzeigerPredict()'s own design) - so the training matrix and
times are carried in the returned object as well as the SPC result.
fitZeitZeiger( trainMat, trainTime, period = 24, nKnots = 3, nTime = 10, sumabsv = 1, orth = TRUE, useSpc = TRUE )fitZeitZeiger( trainMat, trainTime, period = 24, nKnots = 3, nTime = 10, sumabsv = 1, orth = TRUE, useSpc = TRUE )
trainMat |
Numeric matrix of expression, genes x samples. |
trainTime |
Numeric vector of known circadian phase in |
period |
Length of the cycle (default 24, for hours). |
nKnots |
Number of internal knots for the periodic smoothing spline, passed
to |
nTime |
Number of time-points used to discretise each gene's fitted curve
before computing sparse principal components, passed to |
sumabsv |
L1-constraint on the sparse principal components, passed to
|
orth |
Logical; require the components be orthogonal, passed to
|
useSpc |
Logical; use |
An object of class "zeitZeigerFit": a list with the training gene names
(genes), the transposed training matrix (xTrain) and fractional training
time (timeTrain) that zeitzeigerPredict() needs at prediction time, the SPC
result (spcResult), nKnots, and period.
if (requireNamespace("zeitzeiger", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitZeitZeiger(trainMat, time) length(fit$genes) }if (requireNamespace("zeitzeiger", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitZeitZeiger(trainMat, time) length(fit$genes) }
Wraps babelgene::orthologs() (AGENTS.md Section 4 Step 1.2), restricted to
unambiguous one-to-one orthologs: a gene with zero or more than one candidate
ortholog is dropped (ortholog = NA, status records why), never guessed at
via a many-to-many mapping.
mapOrthologs(genes, fromSpecies, toSpecies)mapOrthologs(genes, fromSpecies, toSpecies)
genes |
Character vector of gene symbols in |
fromSpecies, toSpecies
|
|
babelgene's bundled symbol table is a HomoloGene-vintage snapshot and can lag
current HGNC symbols - verified directly for a gene this package cares about:
the human ortholog of mouse Arntl is stored there under the legacy symbol
"BMAL1", not the current official "ARNTL", so querying "ARNTL" on its own
returns zero rows. To avoid that silently dropping known circadian genes, every
query also tries the gene's known aliases from circaTime's own alias table
(the same one standardizeGeneIds() uses) in a single batched lookup, and
results are relabelled back to the originally requested symbol.
A data.frame, one row per unique input gene (input order preserved):
input_gene, ortholog (NA if none), support_n (number of supporting
ortholog databases, NA if none), status ("ok", "not_found", or
"ambiguous"). A warning() reports how many genes were dropped and why.
mapOrthologs(c("Arntl", "Clock", "Per1"), fromSpecies = "Mmusculus", toSpecies = "Hsapiens")mapOrthologs(c("Arntl", "Clock", "Per1"), fromSpecies = "Mmusculus", toSpecies = "Hsapiens")
The diagnostic for users who are not circadian researchers (AGENTS.md Section 6):
given a design such as ~ condition, is the samples' circadian phase distributed
unevenly across that design's terms? If it is, rhythmic genes will show apparent
condition effects that are really collection-time effects.
phaseConfounding( se, design, phase_col = "phase", flag_col = "phase_flag", method = NULL, tissue = NULL, species = "Mmusculus", nPerm = 1000, period = 24, seed = NULL )phaseConfounding( se, design, phase_col = "phase", flag_col = "phase_flag", method = NULL, tissue = NULL, species = "Mmusculus", nPerm = 1000, period = 24, seed = NULL )
se |
A |
design |
A one-sided formula naming |
phase_col |
Name of the |
flag_col |
Name of the |
method, tissue, species
|
Passed to |
nPerm |
Number of permutations (default 1000). |
period |
Length of the cycle (default 24). |
seed |
Optional integer seed; the global RNG state is left unchanged either way (AGENTS.md Section 10). |
Phase is only a confounder when it is unbalanced across a term's levels. Phase that varies randomly within every group inflates residual variance but does not bias the term's coefficient, so a non-significant result here is genuinely reassuring about bias (though not about power).
A data.frame with one row per tested term: term, type
("categorical"/"numeric"), n (samples used), n_groups (NA for numeric
terms), statistic (Watson-Williams F, or circular-linear correlation),
effect_size (for categorical terms the largest pairwise circular distance
between group mean phases, in the same units as period, i.e. directly
readable as hours; for numeric terms the circular-linear correlation itself),
and p_value.
Two attributes carry the supporting detail: "groupMeans", a data.frame
(term, group, n, circular_mean) of per-group mean phases for the
categorical terms, and "nDropped", the number of samples excluded for having
an NA or non-"ok" phase.
This tests the association between phase and each design term. It does not
test which genes are rhythmic - rhythmicity testing on time-labelled data is a
solved problem and an explicit non-goal (AGENTS.md Section 3); hand the phased
object to limorhyde/dryR for that. It also does not adjust or correct the
design; it reports the imbalance so the analyst can decide whether to add phase
as a covariate, which is a modelling decision, not something to apply silently.
Both tests are circular (AGENTS.md 7.1) and both take their null distribution from permuting phase against the term, which avoids the von Mises and equal-concentration assumptions the parametric versions rely on:
Categorical terms (factor, character, logical): the Watson-Williams between-group statistic in its F form.
Numeric terms: the Mardia/Johnson-Wehrly circular-linear correlation.
P-values use the add-one estimator (1 + #{perm >= observed}) / (nPerm + 1), so
the smallest reportable value is 1 / (nPerm + 1) and a p-value is never
reported as exactly zero.
set.seed(1) # 20 samples: cases collected around ZT6, controls around ZT14 - a design where # collection time tracks the condition almost perfectly. phase <- c(rnorm(10, mean = 6, sd = 1), rnorm(10, mean = 14, sd = 1)) mat <- matrix(rnorm(100 * 20), nrow = 100) rownames(mat) <- paste0("g", seq_len(100)) colnames(mat) <- paste0("s", seq_len(20)) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) SummarizedExperiment::colData(se)$phase <- wrapPhase(phase) SummarizedExperiment::colData(se)$condition <- rep(c("case", "control"), each = 10) phaseConfounding(se, ~ condition, nPerm = 199, seed = 1)set.seed(1) # 20 samples: cases collected around ZT6, controls around ZT14 - a design where # collection time tracks the condition almost perfectly. phase <- c(rnorm(10, mean = 6, sd = 1), rnorm(10, mean = 14, sd = 1)) mat <- matrix(rnorm(100 * 20), nrow = 100) rownames(mat) <- paste0("g", seq_len(100)) colnames(mat) <- paste0("s", seq_len(20)) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) SummarizedExperiment::colData(se)$phase <- wrapPhase(phase) SummarizedExperiment::colData(se)$condition <- rep(c("case", "control"), each = 10) phaseConfounding(se, ~ condition, nPerm = 199, seed = 1)
Bundles everything estimatePhase() needs to predict phase for new samples
against one (species, tissue, method, platform) reference cohort: the engine's
own fit object, a pool of real held-out circular residuals for
phase_lo/phase_hi (AGENTS.md Section 9 Phase 3 - the same residual-pool
mechanism bootstrapPhaseCI() uses, but precomputed once at training time
rather than live-bootstrapped per call, since refitting e.g. tauFisher dozens of
times per user call would be prohibitively slow), and the core-clock PCA basis +
derived threshold for clock_amp/phase_flag (AGENTS.md 7.4). Built by
inst/scripts/train_reference_models.R from GSE54650/GSE54651 and bundled
internally as .referenceModels (see availableReferenceModels() to list what
is covered).
## S4 method for signature 'PhaseReferenceModel' show(object)## S4 method for signature 'PhaseReferenceModel' show(object)
object |
A |
An S4 object of class PhaseReferenceModel.
species,tissue,methodCharacter scalars identifying this model. Only
species = "Mmusculus" is currently bundled - see estimatePhase()'s
documentation for why (AGENTS.md 7.6 gap).
platformCharacter scalar: the data type this model was trained on, e.g.
"microarray" (GSE54650) or "rnaseq" (GSE54651, AGENTS.md Section 5 Step
2.3). Defaults to "microarray" for models built before this slot existed,
so old new("PhaseReferenceModel", ...) calls that omit it keep working.
fitThe engine's own fit object (e.g. from fitMolecularTimetable() or
fitTauFisher()), passed as-is to that engine's predict function. Typed
ANY rather than list: both fit objects are S3-classed lists
("molecularTimetableFit"/"tauFisherFit"), and S4's validObject()
rejects an S3-classed object against a list-typed slot even though
is.list() on it is TRUE - the engines' own predict functions require
that S3 class via inherits(), so it must survive storage in this slot
unchanged.
residualPoolNumeric vector of real held-out circular residuals
(circularDiff(pred, truth)) pooled from tissues OTHER than this model's own
(leave-tissue-out, same principle as benchmarks/run_phase3_calibration.R).
ampCenter,ampScaleNamed numeric vectors: the per-gene z-score
center/scale from clockAmplitude()'s PCA fit on this cohort.
ampRotationNumeric matrix (genes x nPC): the PCA rotation from
clockAmplitude() on this cohort.
ampGenesCharacter vector, the core clock genes (row order matching
ampCenter/ampScale/ampRotation) used for this cohort's amplitude PCA.
ampThresholdNumeric scalar: the clockAmplitude() cutoff below which
phase_flag becomes "low_amplitude" (from deriveAmplitudeThreshold()).
cvMAENumeric scalar: this exact (species, tissue, method, platform)
combination's OWN held-out cross-validated circular MAE, in hours, measured
during training. Distinct from mean(abs(residualPool)), which summarises
other tissues' residuals (the pool is leave-tissue-out) and therefore
cannot say how well this engine does on this particular tissue. That
distinction is the whole point of the slot: estimatePhase()'s
method = "consensus" weights each engine by 1/cvMAE, which is only
tissue-aware because this number is. NA_real_ if not measured.
trainedOnCharacter scalar, a citation/provenance string.
rm <- methods::new("PhaseReferenceModel", species = "Mmusculus", tissue = "Liv", method = "timetable", fit = list(genes = c("Arntl", "Per1")), residualPool = rnorm(50, sd = 1.5), ampCenter = c(Arntl = 0, Per1 = 0), ampScale = c(Arntl = 1, Per1 = 1), ampRotation = matrix(c(0.7, 0.7, 0.7, -0.7), nrow = 2), ampGenes = c("Arntl", "Per1"), ampThreshold = 2, trainedOn = "toy example, not a real fit" ) rmrm <- methods::new("PhaseReferenceModel", species = "Mmusculus", tissue = "Liv", method = "timetable", fit = list(genes = c("Arntl", "Per1")), residualPool = rnorm(50, sd = 1.5), ampCenter = c(Arntl = 0, Per1 = 0), ampScale = c(Arntl = 1, Per1 = 1), ampRotation = matrix(c(0.7, 0.7, 0.7, -0.7), nrow = 2), ampGenes = c("Arntl", "Per1"), ampThreshold = 2, trainedOn = "toy example, not a real fit" ) rm
Base-graphics plot of empirical vs. nominal coverage (AGENTS.md Phase 3 task 3: "Plot it"), with the identity line marking perfect calibration. Points below the line mean the stated interval is too narrow (overconfident); points above mean it is too wide (underconfident, but honest).
plotCalibration(calCurve)plotCalibration(calCurve)
calCurve |
Output of |
calCurve, invisibly. Called for its plotting side effect.
calCurve <- data.frame(nominal = c(0.5, 0.8, 0.95), empirical = c(0.53, 0.81, 0.96)) plotCalibration(calCurve)calCurve <- data.frame(nominal = c(0.5, 0.8, 0.95), empirical = c(0.53, 0.81, 0.96)) plotCalibration(calCurve)
For each test sample, finds the phase (on a fine grid over one cycle) whose per-gene harmonic-fit template best correlates with the sample's observed expression of the fitted model's time-indicating genes.
predictMolecularTimetable(fit, testMat, gridN = 288, minGenes = 5)predictMolecularTimetable(fit, testMat, gridN = 288, minGenes = 5)
fit |
A |
testMat |
Numeric matrix of expression, genes x samples, on the same scale
as the training matrix. Only genes overlapping |
gridN |
Number of grid points spanning one cycle to search over (default 288, i.e. 5-minute resolution over 24h). |
minGenes |
Minimum number of overlapping time-indicating genes required to
attempt a prediction for a sample; below this, the sample gets |
A data.frame with one row per column of testMat: phase (numeric,
[0, period), or NA if prediction failed), and confidence (the maximum
per-sample correlation achieved across the grid, or NA).
set.seed(1) nGenes <- 20; nSamples <- 24 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.1) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitMolecularTimetable(trainMat, trainTime, nGenes = 10, minCor = 0.3) pred <- predictMolecularTimetable(fit, trainMat) head(pred)set.seed(1) nGenes <- 20; nSamples <- 24 trainTime <- runif(nSamples, 0, 24) theta <- trainTime / 24 * 2 * pi truePhase <- runif(nGenes, 0, 2 * pi) trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - truePhase[i]) + rnorm(nSamples, sd = 0.1) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitMolecularTimetable(trainMat, trainTime, nGenes = 10, minCor = 0.3) pred <- predictMolecularTimetable(fit, trainMat) head(pred)
Calls the authors' own tauFisher::test_tauFisher(). Unlike train_tauFisher(),
this step is deterministic given the fitted model (no RNG use).
predictTauFisher(fit, testMat, minGenes = 2)predictTauFisher(fit, testMat, minGenes = 2)
fit |
A |
testMat |
Numeric matrix of expression, genes x samples, same scale as the
training matrix. Prediction requires every one of |
minGenes |
Minimum number of |
A data.frame with one row per column of testMat: sample, phase
(integer hour in [0, 24), from nnet::multinom's discrete class prediction,
or NA if prediction was not attempted), and confidence (NA_real_ for
every row - tauFisher has no native per-sample confidence signal; calibrated
uncertainty is Phase 3 scope, AGENTS.md Section 9).
if (requireNamespace("tauFisher", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 100 + 40 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitTauFisher(trainMat, time, seed = 1) pred <- predictTauFisher(fit, trainMat) head(pred) }if (requireNamespace("tauFisher", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 100 + 40 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitTauFisher(trainMat, time, seed = 1) pred <- predictTauFisher(fit, trainMat) head(pred) }
Calls the authors' own TimeSignatR::predTimeStamp(). cv.glmnet's prediction
matches features by column position, not by name, so - like predictZeitZeiger()
this requires every one of fit$genes to be present in testMat's rownames;
if any is missing, every sample gets NA rather than a prediction computed on a
reordered or shrunken gene set.
predictTimeSignatR(fit, testMat)predictTimeSignatR(fit, testMat)
fit |
A |
testMat |
Numeric matrix of expression, genes x samples, same scale as the training matrix. |
A data.frame with one row per column of testMat: sample, phase
(numeric, [0, period), or NA if testMat was missing any fitted gene), and
confidence (NA_real_ for every row - predTimeStamp() returns a point
estimate with no native per-sample confidence signal; calibrated uncertainty is
Phase 3 scope, AGENTS.md Section 9).
if (requireNamespace("TimeSignatR", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) colnames(trainMat) <- paste0("s", seq_len(length(time))) fit <- fitTimeSignatR(trainMat, time, seed = 1) pred <- predictTimeSignatR(fit, trainMat) head(pred) }if (requireNamespace("TimeSignatR", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) colnames(trainMat) <- paste0("s", seq_len(length(time))) fit <- fitTimeSignatR(trainMat, time, seed = 1) pred <- predictTimeSignatR(fit, trainMat) head(pred) }
Calls the authors' own zeitzeiger::zeitzeigerPredict(). zeitzeiger matches
training and test features by column position, not by name, so - unlike
predictMolecularTimetable(), which silently uses whatever fitted genes overlap
testMat - this requires every one of fit$genes to be present in testMat's
rownames; if any is missing, every sample gets NA rather than a prediction
computed on a reordered or shrunken gene set (same all-or-nothing rule as
predictTauFisher()).
predictZeitZeiger(fit, testMat, nSpc = 2, timeRange = seq(0, 1 - 0.01, 0.01))predictZeitZeiger(fit, testMat, nSpc = 2, timeRange = seq(0, 1 - 0.01, 0.01))
fit |
A |
testMat |
Numeric matrix of expression, genes x samples, same scale as the training matrix. |
nSpc |
Number of sparse principal components to use for prediction, passed
to |
timeRange |
Vector of fractional-time grid points searched for the
maximum-likelihood initial value, passed through to |
A data.frame with one row per column of testMat: sample, phase
(numeric, [0, period), or NA if testMat was missing any fitted gene), and
confidence (NA_real_ for every row - zeitzeiger's per-sample likelihood
surface is not converted to a calibrated confidence score here; calibrated
uncertainty is Phase 3 scope, AGENTS.md Section 9).
if (requireNamespace("zeitzeiger", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitZeitZeiger(trainMat, time) pred <- predictZeitZeiger(fit, trainMat) head(pred) }if (requireNamespace("zeitzeiger", quietly = TRUE)) { set.seed(1) nGenes <- 40; time <- seq(0, 23, by = 1) theta <- time / 24 * 2 * pi trainMat <- t(sapply(seq_len(nGenes), function(i) { 5 + 2 * cos(theta - runif(1, 0, 2 * pi)) + rnorm(length(time), sd = 0.3) })) rownames(trainMat) <- paste0("gene", seq_len(nGenes)) fit <- fitZeitZeiger(trainMat, time) pred <- predictZeitZeiger(fit, trainMat) head(pred) }
Auto-detects the identifier type of rownames(se) per gene (Ensembl, Entrez,
RefSeq, or already gene symbol - AGENTS.md Section 4 Step 1.1) and maps
everything to gene symbol, the internal identifier convention every other
circaTime function uses (estimatePhase(), phaseConfounding(),
clockAmplitude()). Runs in two layers: a dependency-free pass first (Ensembl
version-suffix stripping and the hand-curated circadian clock-gene alias table
also used internally by estimatePhase()), then, only for whatever remains
unresolved, AnnotationDbi::mapIds() against org.Mm.eg.db/org.Hs.eg.db if
installed (Suggests, not required for the common already-symbol case - see
notes/decisions.md Section 8).
standardizeGeneIds(se, species = "Mmusculus", map_aliases = TRUE)standardizeGeneIds(se, species = "Mmusculus", map_aliases = TRUE)
se |
A |
species |
|
map_aliases |
Logical, default |
This does not collapse rows that map to the same symbol (e.g. two RefSeq
transcripts of one gene) - rownames(se) may contain duplicates afterward,
same as an unstandardized SummarizedExperiment can already have. Functions
that need a unique gene set (e.g. estimatePhase()'s internal
.prepareExpressionForModel()) already collapse duplicates by row mean at
that point.
se with rownames(se) replaced by standardized gene symbols.
metadata(se)$circaTime_id_mapping is set to a one-row-per-input-gene
data.frame: input_id, id_type ("symbol"/"ensembl"/"entrez"/
"refseq"), mapped_symbol, resolved_by ("unchanged", "free_layer"
– Ensembl version-stripping, the hardcoded core-clock-gene Ensembl table, or
the alias table –, "annotationdbi", or "unresolved"). A summary is also
emitted via message(): how many genes were resolved, and by which layer.
mat <- matrix(rnorm(9), nrow = 3, dimnames = list(c("Arntl", "Bmal1", "Clock"), paste0("s", 1:3))) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) se <- standardizeGeneIds(se, species = "Mmusculus") rownames(se)mat <- matrix(rnorm(9), nrow = 3, dimnames = list(c("Arntl", "Bmal1", "Clock"), paste0("s", 1:3))) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) se <- standardizeGeneIds(se, species = "Mmusculus") rownames(se)
Fits each engine on one SummarizedExperiment and scores it on a different
one (AGENTS.md Section 5 Step 2.2), e.g. training on microarray data and
predicting on RNA-seq data from the same tissue/study. This is a distinct
question from benchmarkPhase()'s cross-validation (how well an engine
predicts held-out samples from the same dataset it was trained on) and
benchmarkPhase() has no way to ask it - there is exactly one train/test split
here, fixed by which dataset a sample belongs to, not nFolds random ones. This
function is not a parallel benchmarking harness, though: it reuses the same
engine registry, gene-name-based fit/predict matching (each engine's own
predict function already intersects fit$genes against rownames(testMat),
so trainSe/testSe do not need to share an identical gene panel), and
circular-error metrics that benchmarkPhase() uses internally.
transferPhase( trainSe, testSe, truth_col, methods = "timetable", train_assay_name = NULL, test_assay_name = NULL, period = 24, engineArgs = list() )transferPhase( trainSe, testSe, truth_col, methods = "timetable", train_assay_name = NULL, test_assay_name = NULL, period = 24, engineArgs = list() )
trainSe, testSe
|
|
truth_col |
Name of the |
methods |
Character vector of engine names; see |
train_assay_name, test_assay_name
|
Assay to use in |
period |
Length of the cycle (default 24). Passed through to every engine
whose |
engineArgs |
Named list, keyed by method, of extra arguments to pass to
that engine's |
A list with two elements: summary, a data.frame with one row per
method (method, circular_mae, circular_cor, n, n_failed); and
residuals, a data.frame with one row per (method, testSe sample)
(method, sample, truth, pred, circular_error, aligned, failed).
If an engine's fit/predict call actually errors, that method's
predictions are all NA and n_failed records how many test samples were
lost, exactly as benchmarkPhase() does for a failed fold. Separately, three
of the four wrapped engines (predictZeitZeiger(), predictTimeSignatR(),
and predictTauFisher()) require the entire fitted gene set to be present
in testMat and silently return all-NA phase (no error at all) when it is
not - only predictMolecularTimetable() tolerates partial gene-panel overlap
(see their own docs). That case shows up here as n = 0 with n_failed = 0
(not counted as a "failure" - nothing errored), not n_failed = n; check n
as well as n_failed when auditing a transfer result, and see
notes/decisions.md for a worked example against real cross-platform data.
This never fabricates a result for a failed sample (AGENTS.md Section 0.3).
set.seed(1) makeSe <- function(n) { time <- seq(0, 23, length.out = n) theta <- time / 24 * 2 * pi mat <- rbind( g1 = 5 + 2 * cos(theta) + rnorm(n, sd = 0.2), g2 = 5 + 2 * sin(theta) + rnorm(n, sd = 0.2) ) colnames(mat) <- paste0("s", seq_len(n)) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) SummarizedExperiment::colData(se)$ZT <- time se } tb <- transferPhase(makeSe(12), makeSe(8), truth_col = "ZT", methods = "timetable") tb$summaryset.seed(1) makeSe <- function(n) { time <- seq(0, 23, length.out = n) theta <- time / 24 * 2 * pi mat <- rbind( g1 = 5 + 2 * cos(theta) + rnorm(n, sd = 0.2), g2 = 5 + 2 * sin(theta) + rnorm(n, sd = 0.2) ) colnames(mat) <- paste0("s", seq_len(n)) se <- SummarizedExperiment::SummarizedExperiment(assays = list(logExpr = mat)) SummarizedExperiment::colData(se)$ZT <- time se } tb <- transferPhase(makeSe(12), makeSe(8), truth_col = "ZT", methods = "timetable") tb$summary
[0, period)
Wrap a phase value onto [0, period)
wrapPhase(x, period = 24)wrapPhase(x, period = 24)
x |
Numeric vector of phase values, in the same units as |
period |
Length of the cycle (default 24, for hours). |
Numeric vector, each element wrapped into [0, period).
wrapPhase(c(-1, 0, 23.5, 24, 25.5))wrapPhase(c(-1, 0, 23.5, 24, 25.5))