Tens of thousands of public RNA-seq samples have no recorded
collection time, yet circadian phase — the time of sample collection
relative to the internal clock — can often be recovered computationally
from expression data alone. Several published methods do this, but each
lives in its own repository, takes its own input format, and reports its
own notion of confidence (usually none at all). circaTime
wraps published phase-inference methods behind a single
SummarizedExperiment-based interface, adds a benchmarking
harness so methods can be compared on the same labelled ground truth,
and adds a calibrated uncertainty layer — a confidence interval and a
quality flag — rather than a bare point estimate.
This vignette walks through the whole pipeline: a raw
SummarizedExperiment, the call that adds phase estimates to
it, how to check those estimates against known ground truth, what the
added confidence columns mean, and how a phased
SummarizedExperiment hands off to downstream rhythmicity
analysis.
The example data bundled with the package is a real subset of GSE54650 (Zhang et al. 2014, PNAS): mouse liver, 24 samples spanning one full circadian cycle at 2-hour resolution, full gene set. Real published data was used here specifically so every number in this vignette is a real result, not a simulation.
se <- readRDS(system.file("extdata", "example_se.rds", package = "circaTime"))
se
#> class: SummarizedExperiment
#> dim: 22105 24
#> metadata(4): source platform mapping_rate note
#> assays(2): expr logExpr
#> rownames(22105): 0610007P14Rik 0610009B22Rik ... Zzef1 Zzz3
#> rowData names(0):
#> colnames(24): GSM1321182 GSM1321183 ... GSM1321204 GSM1321205
#> colData names(6): sample_title geo_accession ... phase donorcolData(se) happens to carry a phase column
here because GSE54650’s collection time is known — that is what lets
this vignette also demonstrate benchmarking below. Most real-world input
to estimatePhase() will not have this column; that
is precisely the value the package provides.
head(colData(se))
#> DataFrame with 6 rows and 6 columns
#> sample_title geo_accession tissue ct phase
#> <character> <character> <character> <numeric> <numeric>
#> GSM1321182 Liv_CT18 GSM1321182 Liv 18 18
#> GSM1321183 Liv_CT20 GSM1321183 Liv 20 20
#> GSM1321184 Liv_CT22 GSM1321184 Liv 22 22
#> GSM1321185 Liv_CT24 GSM1321185 Liv 24 0
#> GSM1321186 Liv_CT26 GSM1321186 Liv 26 2
#> GSM1321187 Liv_CT28 GSM1321187 Liv 28 4
#> donor
#> <character>
#> GSM1321182 GSM1321182
#> GSM1321183 GSM1321183
#> GSM1321184 GSM1321184
#> GSM1321185 GSM1321185
#> GSM1321186 GSM1321186
#> GSM1321187 GSM1321187estimatePhase()estimatePhase() is the package’s main entry point. It
requires a tissue (there is no default — applying a model
trained on one tissue to another silently is a bug, not a convenience)
and looks up a matching pre-trained reference model.
availableReferenceModels() lists what is currently
bundled:
availableReferenceModels()
#> species tissue method platform cv_mae n_amp_genes
#> 1 Mmusculus Adr timetable microarray 0.8194444 16
#> 2 Mmusculus Adr taufisher microarray 0.5714286 16
#> 3 Mmusculus Adr consensus_accuracy microarray 0.7109746 16
#> 4 Mmusculus Adr consensus_equal microarray 0.7274306 16
#> 5 Mmusculus Aor timetable microarray 0.9375000 16
#> 6 Mmusculus Aor taufisher microarray 0.4285714 16
#> 7 Mmusculus Aor consensus_accuracy microarray 0.6386457 16
#> 8 Mmusculus Aor consensus_equal microarray 0.7222222 16
#> 9 Mmusculus Bstm timetable microarray 0.6770833 16
#> 10 Mmusculus Bstm taufisher microarray 1.7857143 16
#> 11 Mmusculus Bstm consensus_accuracy microarray 0.7846040 16
#> 12 Mmusculus Bstm consensus_equal microarray 0.9322917 16
#> 13 Mmusculus BFat timetable microarray 0.8159722 16
#> 14 Mmusculus BFat taufisher microarray 0.6428571 16
#> 15 Mmusculus BFat consensus_accuracy microarray 0.7355824 16
#> 16 Mmusculus BFat consensus_equal microarray 0.7447917 16
#> 17 Mmusculus Cer timetable microarray 0.9618056 16
#> 18 Mmusculus Cer taufisher microarray 1.2857143 16
#> 19 Mmusculus Cer consensus_accuracy microarray 1.0383770 16
#> 20 Mmusculus Cer consensus_equal microarray 1.0520833 16
#> 21 Mmusculus Hrt timetable microarray 0.7673611 16
#> 22 Mmusculus Hrt taufisher microarray 0.6428571 16
#> 23 Mmusculus Hrt consensus_accuracy microarray 0.6759288 16
#> 24 Mmusculus Hrt consensus_equal microarray 0.6822917 16
#> 25 Mmusculus Hyp timetable microarray 1.0138889 16
#> 26 Mmusculus Hyp taufisher microarray 1.2142857 16
#> 27 Mmusculus Hyp consensus_accuracy microarray 1.0610403 16
#> 28 Mmusculus Hyp consensus_equal microarray 1.0694444 16
#> 29 Mmusculus Kid timetable microarray 0.6527778 16
#> 30 Mmusculus Kid taufisher microarray 0.5714286 16
#> 31 Mmusculus Kid consensus_accuracy microarray 0.6259304 16
#> 32 Mmusculus Kid consensus_equal microarray 0.6267361 16
#> 33 Mmusculus Liv timetable microarray 0.7500000 16
#> 34 Mmusculus Liv taufisher microarray 0.7142857 16
#> 35 Mmusculus Liv consensus_accuracy microarray 0.6966251 16
#> 36 Mmusculus Liv consensus_equal microarray 0.6979167 16
#> 37 Mmusculus Lun timetable microarray 0.4895833 16
#> 38 Mmusculus Lun taufisher microarray 0.5000000 16
#> 39 Mmusculus Lun consensus_accuracy microarray 0.4655359 16
#> 40 Mmusculus Lun consensus_equal microarray 0.4652778 16
#> 41 Mmusculus Mus timetable microarray 0.8854167 16
#> 42 Mmusculus Mus taufisher microarray 0.5714286 16
#> 43 Mmusculus Mus consensus_accuracy microarray 0.8093408 16
#> 44 Mmusculus Mus consensus_equal microarray 0.8194444 16
#> 45 Mmusculus WFat timetable microarray 0.6076389 16
#> 46 Mmusculus WFat taufisher microarray 0.6428571 16
#> 47 Mmusculus WFat consensus_accuracy microarray 0.5378838 16
#> 48 Mmusculus WFat consensus_equal microarray 0.5399306 16
#> 49 Mmusculus Adr timetable rnaseq 5.1875000 16
#> 50 Mmusculus Aor timetable rnaseq 3.7604167 16
#> 51 Mmusculus BFat timetable rnaseq 4.5000000 16
#> 52 Mmusculus Bstm timetable rnaseq 2.8854167 16
#> 53 Mmusculus Cer timetable rnaseq 6.1666667 16
#> 54 Mmusculus Hrt timetable rnaseq 3.2812500 16
#> 55 Mmusculus Hyp timetable rnaseq 4.7187500 16
#> 56 Mmusculus Kid timetable rnaseq 1.1666667 16
#> 57 Mmusculus Liv timetable rnaseq 3.9583333 16
#> 58 Mmusculus Lun timetable rnaseq 1.8333333 16
#> 59 Mmusculus Mus timetable rnaseq 3.2604167 16
#> 60 Mmusculus WFat timetable rnaseq 4.3958333 16
#> n_residual_draws
#> 1 264
#> 2 154
#> 3 264
#> 4 264
#> 5 264
#> 6 154
#> 7 264
#> 8 264
#> 9 264
#> 10 154
#> 11 264
#> 12 264
#> 13 264
#> 14 154
#> 15 264
#> 16 264
#> 17 264
#> 18 154
#> 19 264
#> 20 264
#> 21 264
#> 22 154
#> 23 264
#> 24 264
#> 25 264
#> 26 154
#> 27 264
#> 28 264
#> 29 264
#> 30 154
#> 31 264
#> 32 264
#> 33 264
#> 34 154
#> 35 264
#> 36 264
#> 37 264
#> 38 154
#> 39 264
#> 40 264
#> 41 264
#> 42 154
#> 43 264
#> 44 264
#> 45 264
#> 46 154
#> 47 264
#> 48 264
#> 49 88
#> 50 88
#> 51 88
#> 52 88
#> 53 88
#> 54 88
#> 55 88
#> 56 88
#> 57 88
#> 58 88
#> 59 88
#> 60 88phased <- estimatePhase(se, method = "timetable", tissue = "Liv")
colData(phased)[, c("phase", "phase_lo", "phase_hi", "clock_amp", "phase_flag")]
#> DataFrame with 24 rows and 5 columns
#> phase phase_lo phase_hi clock_amp phase_flag
#> <numeric> <numeric> <numeric> <numeric> <factor>
#> GSM1321182 17.66667 16.5250 18.833333 2.18624 low_amplitude
#> GSM1321183 20.16667 19.0250 21.333333 2.99232 ok
#> GSM1321184 22.33333 21.1917 23.500000 4.11958 ok
#> GSM1321185 23.50000 22.3583 0.666667 5.30768 ok
#> GSM1321186 1.91667 0.7750 3.083333 5.15771 ok
#> ... ... ... ... ... ...
#> GSM1321201 8.83333 7.69167 10.0000 4.62902 ok
#> GSM1321202 11.08333 9.94167 12.2500 3.89357 ok
#> GSM1321203 12.16667 11.02500 13.3333 3.40966 ok
#> GSM1321204 15.00000 13.85833 16.1667 2.11778 low_amplitude
#> GSM1321205 16.75000 15.60833 17.9167 2.28560 low_amplitudeFive columns are added: phase (the point estimate, hours
in [0, 24)), phase_lo/phase_hi (a
confidence interval — see below), clock_amp (core-clock
amplitude — see below), and phase_flag (a quality flag:
"ok", "low_amplitude",
"gene_coverage", or "failed").
estimatePhase() never hard-errors on a difficult individual
sample; a sample that cannot be scored gets phase = NA and
an informative flag instead. Only malformed input (an unrecognised
tissue/species/method combination, or an assay with no gene identifiers)
is a hard error.
method = "consensus"No engine wins everywhere — across the package’s full 12-tissue
GSE54650 benchmark (benchmarks/run_consensus_benchmark.R),
three different engines each win in at least three tissues. Rather than
forcing that choice on you, method = "consensus" runs every
engine bundled for that tissue (for mouse microarray: Molecular
Timetable and tauFisher) and combines their estimates on the circle:
The two chunks below are shown but not executed here:
method = "consensus" combines timetable with
tauFisher, and tauFisher is a GitHub-only
package (Suggests, not installed by Bioconductor’s or
r-universe’s build system). Install it locally
(remotes::install_github("micnngo/tauFisher")) and re-run
these chunks to see the real numbers.
consensus <- estimatePhase(se, method = "consensus", tissue = "Liv")
colData(consensus)[, c("phase", "phase_lo", "phase_hi", "phase_flag")]The combination is not a black box. Everything that went into it is returned:
report <- S4Vectors::metadata(consensus)$circaTime_consensus
report$stored_weights # the trained prior, per engine
report$per_engine_phase # what each engine said, per sample
round(report$concentration, 3) # how strongly the engines agreed, in [0, 1]Two weighting schemes are available. The default,
consensus_weights = "accuracy", weights each engine by
1/cvMAE using that engine’s own held-out error on
this specific tissue (visible as cv_mae in
availableReferenceModels()), so an engine that does well on
liver but poorly on brainstem is weighted differently in each.
consensus_weights = "equal" is a plain circular mean,
useful as the honest baseline for whether the accuracy weighting is
buying anything. These are two separately trained models, not one model
with a switch — each carries its own confidence interval, measured under
its own weighting.
If one engine’s genes are missing from your data, its weight drops to
zero and the rest are renormalised, rather than it contributing an
unusable estimate; the result is only flagged
"gene_coverage" if no engine had coverage.
Benchmarked with real donor-split cross-validation across all 12
GSE54650 tissues (benchmarks/run_consensus_benchmark.R):
pooled mean circular MAE was 0.723h (accuracy-weighted consensus),
0.726h (equal-weighted), 0.738h (Molecular Timetable alone), 0.804h
(tauFisher alone) — both consensus schemes beat committing to either
single engine everywhere, though the two weighting schemes are close
enough to each other (0.723h vs. 0.726h) that this is not a strong
argument for tissue-aware weighting specifically over the simpler
equal-weighted baseline. See notes/decisions.md Section 14
for the full accounting, including the oracle
(best-single-engine-per-tissue, not achievable in practice)
comparison.
Real input data does not always arrive as mouse gene symbols.
standardizeGeneIds() auto-detects whether
rownames(se) are gene symbols, Ensembl IDs, Entrez IDs, or
RefSeq IDs (per gene, so mixed input is fine), and maps everything to
gene symbol — the identifier convention every other function in this
package uses:
mixedIds <- SummarizedExperiment(
assays = list(logExpr = matrix(rnorm(9), nrow = 3,
dimnames = list(c("Arntl", "ENSMUSG00000029238", "11461"), c("s1", "s2", "s3"))))
)
mixedIds <- standardizeGeneIds(mixedIds, species = "Mmusculus")
#> 'select()' returned 1:1 mapping between keys and columns
#> standardizeGeneIds(): 3 genes; 1 resolved by the free layer, 1 by AnnotationDbi, 0 unresolved (3/3 = 100.0% mapped).
rownames(mixedIds) # Arntl (already a symbol), Clock (Ensembl), Actb (Entrez)
#> [1] "Arntl" "Clock" "Actb"Two layers do the work. A dependency-free layer runs first: it strips
Ensembl version suffixes and resolves a small, hand-curated table of
common circadian clock-gene synonyms
(e.g. BMAL1/MOP3 → ARNTL,
DEC1 → BHLHE40, CHRONO/
C1orf51 → CIART) — curated by hand from
HGNC/MGI “previous symbol” and “alias” fields for the core clock genes
this package’s bundled models and QC panel actually use, not pulled
programmatically from a live database (so it is small, auditable in
R/gene-mapping.R, and does not need a network call).
Anything the free layer cannot resolve — a real Ensembl/Entrez/RefSeq
ID, or a symbol outside that curated set — falls through to
AnnotationDbi::mapIds() against
org.Mm.eg.db/org.Hs.eg.db if installed
(Suggests, not required for the common already-symbol
case).
Cross-species mapping is a separate, harder problem from
cross-ID-type mapping within one species, and
mapOrthologs() handles it via the babelgene
package’s curated one-to-one mouse↔︎human ortholog tables — a gene with
no unambiguous ortholog is dropped (NA), not guessed
at:
mapOrthologs(c("Arntl", "Clock", "Per1"), fromSpecies = "Mmusculus", toSpecies = "Hsapiens")
#> input_gene ortholog support_n status
#> 1 Arntl ARNTL 11 ok
#> 2 Clock CLOCK 12 ok
#> 3 Per1 PER1 11 okestimatePhase() wires this in directly via
data_species: set it when se‘s genes are in a
different species’ convention than the bundled model (currently
mouse-only) you are calling against, e.g. human input data against the
bundled mouse model. map_orthologs = TRUE is the default
whenever the two differ, and a mapping summary is recorded at
metadata(se)$circaTime_ortholog_mapping, the same way
circaTime_gene_mapping already records the within-species
alias/ID-type report.
Because this example happens to carry real ground truth, we can check
point-estimate accuracy directly with benchmarkPhase(),
which cross-validates a method with held-out folds split by donor (never
by sample, so no individual’s data leaks across train/test):
bm <- benchmarkPhase(se, truth_col = "phase", methods = "timetable", nFolds = 5, seed = 1)
bm$summary
#> method circular_mae circular_cor n n_failed
#> 1 timetable 0.75 -0.5898503 24 0circular_mae is mean absolute error on a 24-hour circle
— computed as pmin(abs(d), 24 - abs(d)), since 23.5h and
0.5h are one hour apart, not 23. Ordinary (non-circular) statistics are
not meaningful on phase and are never used anywhere in this package.
This result matches what was reported for liver in the package’s own
benchmark suite (circular_mae \(\approx\) 0.98h across the full 12-tissue
GSE54650 atlas; see the package’s benchmarks/results/ for
the complete, reproducible table).
benchmarkPhase() answers “how accurate is this engine on
held-out samples from this dataset?” A different, complementary
question is transfer: “if I train on one dataset/platform, how accurate
is it on a different one?” — e.g. training on microarray data
and scoring RNA-seq samples of the same tissue.
transferPhase() answers that, reusing the same engine
registry and circular-error metrics as benchmarkPhase(),
just with the train/test split fixed by which dataset a sample came from
instead of drawn as cross-validation folds:
seMicroarray <- circaTimeData::circaTimeData("GSE54650") # same study, microarray
seRnaseq <- circaTimeData::circaTimeData("GSE54651") # same study, RNA-seq
liverMicro <- seMicroarray[, SummarizedExperiment::colData(seMicroarray)$tissue == "Liv"]
liverRnaseq <- seRnaseq[, SummarizedExperiment::colData(seRnaseq)$tissue == "Liv"]
tb <- transferPhase(liverMicro, liverRnaseq, truth_col = "phase", methods = "timetable")
tb$summaryphase_lo/phase_hi are not a bootstrap refit
on every call (that would be far too slow for some engines) — they are
built from a pool of real, held-out prediction residuals collected once
when each reference model was trained, added directly to the point
estimate. This was checked directly against real held-out data: a stated
80% interval that used only training-resampling variance (no residual
pool) covered the true phase only about 48% of the time — badly
overconfident. Adding the residual-pool mechanism brought empirical
coverage in line with (in fact slightly above) the nominal level at
every level tested. An uncalibrated interval is worse than none, so this
check was treated as a hard requirement, not an afterthought.
phase_credible_lo/phase_credible_hi is a
second, complementary interval: a Bayesian (von Mises) credible interval
built from the same residual pool as
phase_lo/phase_hi, via
bayesianPhasePosterior(). Where the empirical interval is
the raw quantile spread of the observed residuals, the Bayesian one fits
a smooth parametric distribution to the pool’s concentration and returns
its analytic interval - useful when a caller needs more than interval
endpoints (e.g. bayesianPhasePosterior() also returns a
full posterior density, enabling statements like “what is the
probability the true phase is within 2h of the point estimate”). The two
intervals are expected to be close, not identical: on real held-out
GSE54650 predictions, pooled empirical coverage at nominal 80% was 88.9%
(Bayesian) vs. 87.2% (bootstrap), with mean widths 2.72h vs. 2.78h
(benchmarks/run_bayesian_ci_benchmark.R).
kappa/concentration are recorded in
metadata(se)$circaTime_bayesian.
metadata(phased)$circaTime_bayesian
#> $kappa
#> [1] 13.17283
#>
#> $concentration
#> [1] 0.9612638
#>
#> $level
#> [1] 0.8
colData(phased)[, c("phase", "phase_lo", "phase_hi",
"phase_credible_lo", "phase_credible_hi")]
#> DataFrame with 24 rows and 5 columns
#> phase phase_lo phase_hi phase_credible_lo phase_credible_hi
#> <numeric> <numeric> <numeric> <numeric> <numeric>
#> GSM1321182 17.66667 16.5250 18.833333 16.290974 19.030353
#> GSM1321183 20.16667 19.0250 21.333333 18.790974 21.530353
#> GSM1321184 22.33333 21.1917 23.500000 20.957641 23.697020
#> GSM1321185 23.50000 22.3583 0.666667 22.124308 0.863686
#> GSM1321186 1.91667 0.7750 3.083333 0.540974 3.280353
#> ... ... ... ... ... ...
#> GSM1321201 8.83333 7.69167 10.0000 7.45764 10.1970
#> GSM1321202 11.08333 9.94167 12.2500 9.70764 12.4470
#> GSM1321203 12.16667 11.02500 13.3333 10.79097 13.5304
#> GSM1321204 15.00000 13.85833 16.1667 13.62431 16.3637
#> GSM1321205 16.75000 15.60833 17.9167 15.37431 18.1137clock_amp (phase_flag == "low_amplitude"
when it falls below threshold) has two implementations, selected with
qc_method:
qc_method = "legacy" (default): the Euclidean norm of a
sample’s first two principal components over a fixed
16-gene core clock panel, projected onto the reference cohort’s own PCA
basis (works even for a single-sample call). Honest
caveat: on the two real labelled datasets it was tested against
(the healthy-tissue GSE54650 atlas used above, and a young-vs-aged mouse
stem cell dataset with real circadian-time ground truth), low-amplitude
samples did not show a statistically significant error difference from
high-amplitude ones. Core clock gene expression is reported in
the aging literature to stay close to unaltered while
downstream oscillators desynchronise — so a fixed core-clock
panel may simply be looking at the wrong genes.qc_method = "dynamic" (clockCoherence()):
selects the informative gene panel from the input dataset
itself — each gene’s cosinor fit is tested against a
phase-permutation null, panel = genes clearing a BH-FDR cutoff — rather
than trusting one fixed panel across every dataset. This needs a real
cohort in the call (qc_min_samples, default 10) since the
panel is fit from it; below that it silently falls back to
"legacy"
(metadata(se)$circaTime_qc$method_used records which ran).
On the same young-vs-aged dataset where the legacy flag scored p=0.20,
the dynamic score reaches a selection-corrected p=0.0185 — see
benchmarks/run_dynamic_qc_benchmark.R and
notes/decisions.md Section 15 for the full validation and
why coherence is a self-consistency score (it is scored against
estimatePhase()’s own phase estimate, not independent
ground truth) rather than a general-purpose sample-quality metric.se_dyn <- estimatePhase(se, method = "timetable", tissue = "Liv",
qc_method = "dynamic")
table(colData(se_dyn)$phase_flag)
#>
#> ok low_amplitude gene_coverage failed
#> 24 0 0 0
metadata(se_dyn)$circaTime_qc[c("method_used", "threshold", "panelSize")]
#> $method_used
#> [1] "dynamic"
#>
#> $threshold
#> [1] 0
#>
#> $panelSize
#> all
#> 5860A phased SummarizedExperiment is meant to be handed
straight to existing rhythmicity-detection tools rather than
re-implementing what they already do. The pattern below (illustrative —
limorhyde and dryR are not dependencies of
circaTime and are not installed or executed here) shows the
shape of that handoff: circaTime supplies
colData(se)$phase as a covariate, the downstream tool takes
it from there.
Most users of this package are not circadian researchers. They have a differential expression design — cases versus controls, treated versus untreated — and a reasonable worry: the samples’ collection times were never recorded, so could an unmodelled daily rhythm be driving what looks like a condition effect?
That worry is only justified under a specific condition. Circadian phase biases a contrast when it is unbalanced across the groups being compared. If collection time varies randomly within every group, rhythmic genes inflate residual variance but leave the condition coefficient unbiased. If the cases happened to be collected in the morning and the controls in the afternoon, then every rhythmic gene carries an apparent condition effect that is really a clock effect.
phaseConfounding() tests exactly that: it asks whether
phase is distributed unevenly across each term of a design formula.
# A deliberately confounded design: cases collected around ZT6, controls around ZT14.
set.seed(1)
n <- 20
confPhase <- c(rnorm(10, mean = 6, sd = 1), rnorm(10, mean = 14, sd = 1))
confMat <- matrix(rnorm(100 * n), nrow = 100,
dimnames = list(paste0("g", seq_len(100)), paste0("s", seq_len(n))))
conf <- SummarizedExperiment(assays = list(logExpr = confMat))
colData(conf)$phase <- wrapPhase(confPhase)
colData(conf)$condition <- rep(c("case", "control"), each = 10)
phaseConfounding(conf, ~ condition, nPerm = 999, seed = 1)
#> term type n n_groups statistic effect_size p_value
#> 1 condition categorical 20 2 339.6259 8.132343 0.001The small p-value says the imbalance is real, and
effect_size puts it in interpretable units: the group mean
phases sit roughly eight hours apart. Compare that with a design where
the same samples were collected across the whole day in both arms:
balanced <- conf
colData(balanced)$phase <- rep(seq(0, 22, by = 2), length.out = n)
phaseConfounding(balanced, ~ condition, nPerm = 999, seed = 1)
#> term type n n_groups statistic effect_size p_value
#> 1 condition categorical 20 2 0.577424 4 0.726Here phase is spread evenly through both groups, and the test says so.
Three things are worth being precise about:
effect_size treat them that way. Applying an
ordinary ANOVA or t-test to phase values gets this wrong.If se does not already carry a phase
column, pass method and tissue and
phaseConfounding() will call estimatePhase()
first. Samples whose phase estimate was flagged unreliable
(low_amplitude, gene_coverage,
failed) are excluded from the test rather than allowed to
add noise to it; the number dropped is reported in
attr(result, "nDropped").
This is an early snapshot, and it is worth being explicit about the edges:
species = "Mmusculus") tissues from
GSE54650 are currently bundled; see
availableReferenceModels() for the exact list. A larger,
more tissue/species-diverse reference cohort lives in the separate
companion circaTimeData ExperimentHub package rather than
being bundled directly in circaTime itself.standardizeGeneIds()/mapOrthologs()/estimatePhase(data_species = )
(above) cover ID-type translation and one-to-one mouse↔︎human ortholog
mapping, but there is still no bundled reference model trained
on human data to map onto — an ortholog-mapped human sample is currently
scored against the mouse model, which is unvalidated cross-species
transfer, not a substitute for a real human cohort.phaseConfounding() tests each design term separately.
Interaction terms are skipped, and a full multivariable circular
regression — which would let you ask whether phase is imbalanced across
one term after adjusting for another — is not implemented.sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [3] GenomicRanges_1.65.4 Seqinfo_1.3.2
#> [5] IRanges_2.47.5 S4Vectors_0.51.10
#> [7] BiocGenerics_0.59.12 generics_0.1.4
#> [9] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [11] circaTime_0.99.3 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 SparseArray_1.13.4 RSQLite_3.53.3
#> [4] lattice_0.23-1 magrittr_2.0.5 digest_0.6.39
#> [7] evaluate_1.0.5 grid_4.6.1 fastmap_1.2.0
#> [10] blob_1.3.0 jsonlite_2.0.0 Matrix_1.7-6
#> [13] AnnotationDbi_1.75.2 DBI_1.3.0 BiocManager_1.30.27
#> [16] httr_1.4.9 Biostrings_2.81.9 jquerylib_0.1.4
#> [19] abind_1.4-8 cli_3.6.6 crayon_1.5.3
#> [22] rlang_1.3.0 XVector_0.53.0 bit64_4.8.6
#> [25] withr_3.0.3 cachem_1.1.0 DelayedArray_0.39.8
#> [28] yaml_2.3.12 BiocBaseUtils_1.15.1 otel_0.2.0
#> [31] S4Arrays_1.13.2 tools_4.6.1 memoise_2.0.1
#> [34] dplyr_1.2.1 babelgene_22.9 png_0.1-9
#> [37] buildtools_1.0.0 vctrs_0.7.3 R6_2.6.1
#> [40] lifecycle_1.0.5 KEGGREST_1.53.6 bit_4.6.0
#> [43] pkgconfig_2.0.3 pillar_1.11.1 bslib_0.12.0
#> [46] glue_1.8.1 tidyselect_1.2.1 tibble_3.3.1
#> [49] xfun_0.61 sys_3.4.3 knitr_1.52
#> [52] org.Mm.eg.db_3.23.0 htmltools_0.5.9 rmarkdown_2.32
#> [55] maketools_1.3.2 compiler_4.6.1estimatePhase()method = "consensus"