Movatterモバイル変換

Title:	Changepoint Additive Models for Time Series Omics Data
Version:	0.1.3
Description:	Provides a comprehensive framework for time series omics analysis, integrating changepoint detection, smooth and shape-constrained trends, and uncertainty quantification. It supports gene- and transcript-level inferences, p-value aggregation for improved power, and both case-only and case-control designs. It includes an interactive 'shiny' interface. The methods are described in Yates et al. (2024) <doi:10.1101/2024.12.22.630003>.
License:	GPL (≥ 3)
Encoding:	UTF-8
RoxygenNote:	7.3.2
Imports:	bslib, cli, dplyr, edgeR, ggplot2, grDevices, magrittr,matrixStats, mgcv, mvnfast, parallel, pbmcapply, purrr,RColorBrewer, rlang, scam, shiny, shinyjs, stats, stringr,tidyr, tximport
Depends:	R (≥ 3.5)
URL:	https://l-a-yates.github.io/cpam/,https://github.com/l-a-yates/cpam
Suggests:	knitr, rmarkdown, testthat (≥ 3.0.0)
Config/testthat/edition:	3
BugReports:	https://github.com/l-a-yates/cpam/issues
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2025-03-11 23:41:12 UTC; layates
Author:	Luke Yates [aut, cre, cph], Michael Charleston [aut], Jazmine Humphreys [aut], Steven Smith [aut]
Maintainer:	Luke Yates <luke.yates@utas.edu.au>
Repository:	CRAN
Date/Publication:	2025-03-13 13:20:02 UTC

cpam: Changepoint Additive Models for Time Series Omics Data

Description

Provides a comprehensive framework for time series omics analysis, integrating changepoint detection, smooth and shape-constrained trends, and uncertainty quantification. It supports gene- and transcript-level inferences, p-value aggregation for improved power, and both case-only and case-control designs. It includes an interactive 'shiny' interface. The methods are described in Yates et al. (2024)doi:10.1101/2024.12.22.630003.

Author(s)

Maintainer: Luke Yatesluke.yates@utas.edu.au (ORCID) [copyright holder]

Authors:

• Michael Charleston

• Jazmine Humphreys

• Steven Smith

See Also

Useful links:

• https://l-a-yates.github.io/cpam/

• https://github.com/l-a-yates/cpam

• Report bugs athttps://github.com/l-a-yates/cpam/issues

Pipe operator

Description

Seemagrittr::%>% for details.

Usage

lhs %>% rhs

Arguments

Value

The result of callingrhs(lhs).

Compute normalization factors

Description

Compute normalization factors

Usage

compute_normalization_factors(count_matrix_filtered, normalize)

Arguments

Value

Vector of normalization factors

Compute p-values for each target ID

Description

Compute p-values for each target ID

Usage

compute_p_values(  cpo,  subset = NULL,  p_adj_method = "BH",  gam_method = "REML",  gam_optimizer = "efs",  silent = TRUE)

Arguments

Details

This function computes p-values for each target_id in the supplied cpam object.The p-values are computed from a negative binomial GAM modelwith a thin-plate spline basis function(s) for timeusing themgcv package.

The p-values are stored in the new slotp_table in the cpam object.Ifaggregate_to_gene is set toTRUE (default),the target p-values are aggregated to the gene level using thelancaster method.The columnsp_val_target andp_val_gene store the raw p-values for target- and gene-level, respectively.The function also computes adjusted p-values using thep_adj_method.The default method is "BH" (Benjamini-Hochberg),but any methods supported by the functionp.adjust can be used.The adjusted p-values are stored in the columnsq_val_target andq_val_gene.

Value

an updated cpam object with raw, adjusted, and possibly aggregated p-values stored in the new slot "p_table"

References

Wood, S.N. (2013a) On p-values for smooth components of an extendedgeneralized additive model. Biometrika 100:221-228 doi:10.1093/biomet/ass048

Yi L, Pachter L (2018). aggregation: p-Value Aggregation Methods. R package version 1.0.1,https://CRAN.R-project.org/package=aggregation.

Examples

library(cpam)# load gene-only example cpam objectload(system.file("extdata", "cpo_example.rda", package = "cpam"))# run on a small subset of the example datacpo <- compute_p_values(cpo_example, subset = paste0("g00",1:9))cpo$p_table

Convert count matrix to long format and join with experiment design

Description

Convert count matrix to long format and join with experiment design

Usage

convert_to_long_format(counts_raw, exp_design)

Arguments

Value

Long format data frame

The cpam class

Description

Thecpam class stores data and analysis results for a time seriesomics experiment. This object is generated by theprepare_cpamfunction and contains all necessary data and parameters for downstream analysis.

Details

Acpam object is a list with the following components:

exp_design

A data frame containing the experimental design information, with columns for 'sample', 'time', and potentially other variables.

count_matrix_raw

The original count matrix before filtering.

count_matrix_filtered

The count matrix after filtering low-count genes/transcripts.

target_to_keep

A vector of transcript/gene IDs that passed the filtering criteria.

data_long

A long-format data frame containing all relevant information for each target and sample.

t2g

A transcript-to-gene mapping data frame if provided.

regularize

Logical; whether empirical Bayes regularization of dispersions was used.

overdispersion.prior

Median overdispersion.

model_type

String; the type of design used, either "case-only" or "case-control".

condition_var

String; the column name in exp_design for the condition variable (for case-control models).

case_value

The value of condition_var that indicates the "case" in case-control models.

bootstrap

Logical; whether bootstrap samples (inferential replicates) were used.

nboot

The number of bootstrap samples used, if applicable.

filter

A list containing the filtering function and its arguments used.

gene_level

Logical; whether the analysis was performed at the gene level.

aggregate_to_gene

Logical; whether p-values should be aggregated from transcript to gene level.

times

An ordered vector of unique time points in the experimental design.

num_cores

The number of cores used for parallel computation.

fixed_effects

The formula for fixed effects in the model.

intercept_cc

String; the intercept type for case-control models.

bss

A vector of basis function types used for modelling.

Methods

Objects of classcpam have print and summary methods available.

See Also

prepare_cpam for creating a cpam object.

Examples

# load gene-only example dataload(system.file("extdata", "exp_design_example.rda", package = "cpam"))load(system.file("extdata", "count_matrix_example.rda", package = "cpam"))# Create a cpam object with the example datacpo <- prepare_cpam(exp_design = exp_design_example,                    count_matrix = count_matrix_example,                    gene_level = TRUE)# Print the object structurecpo# Get a summary of the cpam objectsummary(cpo)

Use model selection to estimate changepoints

Description

Use model selection to estimate changepoints

Usage

estimate_changepoint(  cpo,  cps = NULL,  degs_only = TRUE,  deg_threshold = 0.05,  subset = NULL,  sp = NULL,  bss = "tp",  family = c("nb", "gaussian"),  score = "aic",  compute_mvn = TRUE)

Arguments

Details

This function estimates changepoints for each target_id. The assumedtrajectory type for this modelling stage is initially constant followed bya changepoint into thin-plate smoothing spline.

By default, candidate time points are limited to the discrete observed valuesin the series, since, despite the use of smoothing constraints,there is generally insufficient information to infer the timing ofchangepoints beyond the temporal resolution of the data. In any case, thecandidate points can be set manually using thecps argument.

To estimate changepoints, a model is fit for each candidate changepoint andgeneralised cross-validation (GCV, default) or the Akaike InformationCriterion (AIC) are used to select among them. Model-selection uncertaintyis dealt with by computing the one-standard-error rule, which identifies theleast complex model within one standard error of the best scoring model.

Both the minimum and the one-standard-error (default) models are stored in the returnedslot "changepoints" so that either can be used. In addition to these, this function also computes theprobability (denotedp_mvn) that the null model is the best scoring model, using a simulationbased approach based on the multivariate normal model of the pointwisemodel scores.

Given the computational cost of fitting a separate model for each candidatechangepoint, cpam only estimates changepoints for targets associated with'significant' genes at the chosen thresholddeg_threshold.

Value

a cpam object with the estimated changepoint table added to the slot "changepoints"

References

Yates, L. A., S. A. Richards, and B. W. Brook. 2021.Parsimonious model selection using information theory:a modified selection rule.Ecology 102(10):e03475. 10.1002/ecy.3475

Examples

library(cpam)library(dplyr)# load example dataload(system.file("extdata", "exp_design_example.rda", package = "cpam"))load(system.file("extdata", "count_matrix_example.rda", package = "cpam"))cpo <- prepare_cpam(exp_design = exp_design_example,                    count_matrix = count_matrix_example[1:20,],                    gene_level = TRUE,                    num_cores = 1)cpo <- compute_p_values(cpo)cpo <- estimate_changepoint(cpo)cpo$changepoints

Wrapper for estimating dispersions

Description

Wrapper for estimating dispersions

Usage

estimate_dispersions_wrapper(  count_matrix_filtered,  exp_design,  bootstrap,  catch)

Arguments

Value

Estimated dispersions

Handle gene aggregation in long format data

Description

Handle gene aggregation in long format data

Usage

handle_gene_aggregation(data_long, t2g, aggregate_to_gene)

Arguments

Value

Updated long format data with gene ID information

Import count data from files

Description

Import count data from files

Usage

import_count_data(exp_design, t2g, import_type, gene_level, bootstrap)

Arguments

Value

List containing imported data and related variables

Apply the one-standard-error rule for model selection

Description

Apply the one-standard-error rule for model selection

Usage

ose_rule(tab, nse = 1)

Arguments

Details

Implements the one-standard-error rule which selects the most parsimonious modelwithin a specified number of standard errors of the best-scoring model.

Value

The selected model that is most parsimonious within nse standard errors of the best model

Plot clustered targets

Description

Plot clustered targets

Usage

plot_cluster(cpo, res, changepoints, shapes, alpha = 0.1)

Arguments

Details

Plots the fitted trends for a set of targets whose estimated changepoints andshapes are given by the argumentschangepoints andshapes, respectively.

Creates a combined plot showing fitted expression trends for all targets thatshare specified changepoint times and shape patterns. Each line represents onetarget's fitted trajectory, with transparency controlled by alpha.

Value

A ggplot object showing overlaid fitted trends, or NULL if no matchingtargets are found

See Also

results(),plot_cpam()

Examples

library(cpam)# load gene-only example cpam objectload(system.file("extdata", "cpo_example.rda", package = "cpam"))# Generate results tableres_example <- results(cpo_example)# plot all targets with changepoint at timepoint 0 and shape "ilin" (increasing linear)plot_cluster(cpo_example, res_example, changepoints = 2, shapes = "ilin")

Plot fitted changepoint additive models

Description

Plot fitted changepoint additive models

Usage

plot_cpam(  cpo,  gene_id = NULL,  target_id = NULL,  cp_type = c("cp_1se", "cp_min"),  shape_type = "shape1",  bs = "auto",  cp_fix = -1,  facet = FALSE,  sp = NULL,  show_fit = TRUE,  show_data = TRUE,  show_fit_ci = TRUE,  show_data_ci = TRUE,  ci_prob = "se",  remove_null = FALSE,  null_threshold = 0.05,  null_threshold_adj = TRUE,  k_mult = 1.2,  return_fits_only = FALSE,  family = "nb",  common_y_scale = TRUE,  scaled = FALSE,  base_size = 12)

Arguments

Details

Plots the fitted trend and data points for a given gene or target. If a gene IDis supplied, the function will plot all transcripts for that gene.The function can also be used to return the model fit(s) only, which aregamObject objects from themgcv package.

Value

a ggplot object

Examples

library(cpam)# load gene-only example cpam objectload(system.file("extdata", "cpo_example.rda", package = "cpam"))# example geneplot_cpam(cpo_example, gene_id = "g003")# gene with estimated changepoint at timepoint 3plot_cpam(cpo_example, gene_id = "g013")# manually set the changepointplot_cpam(cpo_example, gene_id = "g013", cp_fix = 2)

Prepare a cpam object

Description

Prepare a cpam object

Usage

prepare_cpam(  exp_design,  count_matrix = NULL,  t2g = NULL,  import_type = NULL,  model_type = c("case-only", "case-control"),  bootstrap = TRUE,  filter_fun = "ts_filter",  filter_fun_args = list(min_reads = 5, min_prop = 3/5),  regularize = TRUE,  gene_level = FALSE,  aggregate_to_gene = !gene_level,  condition_var = "condition",  case_value = "treatment",  num_cores = 1,  normalize = TRUE,  fixed_effects = NULL,  intercept_cc = c("1", condition_var))

Arguments

Details

This function prepares a cpam object for analysis.The function loads count data from files or a matrix,filters lowly expressed genes, computes normalisation factors,and estimates dispersions. Many of these steps can be customised or turned off.

When bootstrap samples (inferential replicates) are available, it loads andsummarises these using means, standard errors, and estimated overdispersions.The latter are a measure of quantification uncertainty and they are used torescale the counts which accounts for this uncertainty during the modelling steps.

Value

an object of classcpam-class. The returned object hasmethodsprint andsummary for displaying information.Seecpam-class for details on the structure of the returned object.

References

Pedro L Baldoni, Yunshun Chen, Soroor Hediyeh-zadeh, Yang Liao, Xueyi Dong,Matthew E Ritchie, Wei Shi, Gordon K Smyth,Dividing out quantification uncertainty allows efficient assessment ofdifferential transcript expression with edgeR,Nucleic Acids Research, Volume 52, Issue 3, 9 February 2024,Page e13, https://doi.org/10.1093/nar/gkad1167

Yunshun Chen, Lizhong Chen, Aaron T L Lun, Pedro L Baldoni, Gordon K Smyth,edgeR v4: powerful differential analysis of sequencing data withexpanded functionality and improved support for small countsand larger datasets, Nucleic Acids Research, Volume 53, Issue 2,27 January 2025, https://doi.org/10.1093/nar/gkaf018

Examples

library(cpam)# load gene-only example dataload(system.file("extdata", "exp_design_example.rda", package = "cpam"))load(system.file("extdata", "count_matrix_example.rda", package = "cpam"))cpo <- prepare_cpam(exp_design = exp_design_example,                    count_matrix = count_matrix_example,                    gene_level = TRUE)cpo

Process bootstrap data

Description

Process bootstrap data

Usage

process_bootstrap_data(data_long, catch, boot)

Arguments

Value

Updated data with bootstrap information

Process count matrix when importing is not needed

Description

Process count matrix when importing is not needed

Usage

process_count_matrix(count_matrix, t2g, gene_level)

Arguments

Value

Processed count matrix

Create a results table from a cpam object

Description

Create a results table from a cpam object

Usage

results(  cpo,  p_threshold = 0.05,  p_type = c("p_gam", "p_mvn"),  min_lfc = 0,  min_count = 0,  aggregate_to_gene = cpo$aggregate_to_gene,  add_lfc = TRUE,  add_counts = TRUE,  cp_type = c("cp_1se", "cp_min"),  shape_type = c("shape1", "shape2"),  summarise_to_gene = FALSE,  remove_null_targets = TRUE)

Arguments

Details

This function is usually called aftercompute_p_values(),estimate_changepoint, andselect_shape havebeen run. The function has several useful filters such as adjusted p-valuethresholds, minimum log-fold changes, and minimum counts.

Value

a tibble

Examples

library(cpam)# load gene-only example cpam objectload(system.file("extdata", "cpo_example.rda", package = "cpam"))results(cpo_example)# Add filtersresults(cpo_example, p_threshold = 0.01, min_lfc = 1)

Use model selection to select a shape for each target

Description

Use model selection to select a shape for each target

Usage

select_shape(  cpo,  subset = NULL,  sp = NULL,  bss = c("micv", "mdcx", "cv", "cx", "lin", "tp", "null"),  family = c("nb", "gaussian"),  score = "gcv",  cp_type = c("cp_1se", "cp_min"))

Arguments

Details

The function selects the best shape from a list of candidate shapes for each target.It is typically the last step in the analysis, called after p-values havebeen estimated usingcompute_p_values() and changepoints havebeen estimated usingestimate_changepoint().

Two shape selections are generated. The first selecting among linear,convex and concave shape classes and their monotonic variants (or the shapeset given by bss), and the second selecting among the first options plus an'unconstrained' smooth. The inclusion of the 'unconstrained' type providesthe flexibility to detect targets beyond simpler trends.For computational reasons, as per the changepoint estimation,shapes are selected only for those genes, or their isoforms,identified as significant at the chosen FDR threshold. This isoverridden by providing a subset of target names to thesubset argument,provided these targets have estimated changepoints.

Value

a cpam object with the selected shapes added to the slot "shapes"

Examples

library(cpam)# load example dataload(system.file("extdata", "exp_design_example.rda", package = "cpam"))load(system.file("extdata", "count_matrix_example.rda", package = "cpam"))# Using a small subset of the example datacpo <- prepare_cpam(exp_design = exp_design_example,                    count_matrix = count_matrix_example[1:20,],                    gene_level = TRUE,                    num_cores = 1)cpo <- compute_p_values(cpo)cpo <- estimate_changepoint(cpo)cpo <- select_shape(cpo)cpo$shapes

Simulate p-values using multivariate normal distribution

Description

Simulate p-values using multivariate normal distribution

Usage

simulate_p_mvn(score_table, nsim = 10000, reg = 0.001)

Arguments

Details

Uses Monte Carlo simulation to compute the probability that the null modelis the best scoring model under multivariate normal assumptions.

Value

Probability that the null model is the best scoring model

Summarize bootstrap samples

Description

Internal function to process inferential replicates from tximport objects.

Usage

summarise_bootstraps(txi)

Arguments

Value

A tibble with bootstrap summary statistics

Removes lowly expressed genes

Description

Removes lowly expressed genes

Usage

ts_filter(data, min_reads = 5, min_prop = 3/5)

Arguments

Details

Identifies targets that show strong and consistent expression in at least one timepoint.For each timepoint, the function calculates the proportion of sampleswhere a targets exceedsmin_reads. Targets are retained if they meetthe minimum proportion (min_prop) at any timepoint in the experiment.

Value

a character vector of transcript IDs to keep

Examples

data <- dplyr::tibble(  target_id = rep(paste0("t", 1:3), each = 6),  time = rep(c(0, 4, 8), 6),  counts = c(6,6,6, 0,0,0, 6,0,6, 0,6,6, 6,6,6, 6,0,0))ts_filter(data)

Validate and adjust number of cores

Description

Validate and adjust number of cores

Usage

validate_cores(num_cores)

Arguments

Value

Validated number of cores

Validate input parameters for prepare_cpam

Description

Validate input parameters for prepare_cpam

Usage

validate_inputs(  exp_design,  count_matrix,  t2g,  import_type,  model_type,  condition_var,  case_value,  fixed_effects,  aggregate_to_gene,  gene_level,  import)

Arguments

Value

NULL, throws error if validation fails

Launches a Shiny app to visualise the data and fitted models of a cpam object

Description

Launches a Shiny app to visualise the data and fitted models of a cpam object

Usage

visualise(  cpo,  subset = NULL,  degs_only = TRUE,  deg_threshold = 0.05,  p_type = c("p_gam", "p_mvn"),  shape_type = c("shape1", "shape2"))

Arguments

Value

None (launches Shiny app in browser)

Examples

if (interactive()){# A simple gene-level examplecpo <- cpo_example# Launch visualization with all genesvisualise(cpo, degs_only = FALSE)# Launch with only significant genesvisualise(cpo, deg_threshold = 0.05)# Launch with specific genesvisualise(cpo, subset = c("g001","g002","g003"))}