install.packages("cpam")

Load packages

library(cpam)library(dplyr)library(tidyr)library(stringr)library(ggplot2)

Experimental design

An example experimental design is included in thecpampackage. Since it is case-only design, there are no experimentalconditions beyond time.

# load example dataload(system.file("extdata","exp_design_example.rda",package ="cpam"))exp_design_example

Count matrix

The example count data for are provided as a matrix. Let’s take alook at the first few rows.

# load example dataload(system.file("extdata","count_matrix_example.rda",package ="cpam"))as.data.frame(count_matrix_example)%>% head

Fittingcpam

To fit the models, we first prepare thecpam object,then compute p-values, estimate changepoints, and select the shape foreach gene. In this simple example, simulated data are gene-level, thatis we do not have isoform-level counts. As such, we leave```the transcript-to-gene mapping (t2g) and we setgene_level= T`.

  cpo<-prepare_cpam(exp_design = exp_design_example,count_matrix = count_matrix_example,model_type ="case-only",t2g =NULL,gene_level = T,num_cores =1)# just for the example  cpo<-compute_p_values(cpo)# 6 seconds  cpo<-estimate_changepoint(cpo)# 4 seconds  cpo<-select_shape(cpo)# 5 seconds

We can look at a summary of the fitted cpam object

cpo

If you run the code on your own computer, you can launch the Shinyapp to visualise the results interactively usingvisualise(cpo).

Result tables

The results of the analysis are summarised using theresults function.

results(cpo)

The generated results can be filtered by specifying minimum counts,minimum log-fold changes, and maximum\(p\)-values. For example, to return only thetranscripts with a log-fold change greater than 1, at least 10 counts,and a\(p\)-value less than 0.01, wecan run

results(cpo,min_count =10,min_lfc =1,p_threshold =0.01)

Plotting genes and transcripts

A single gene can be plotted using theplot_cpamfunction. Here we plot the gene g063

plot_cpam(cpo,gene_id ="g063")

The subtitle shows(1,cv) indicating a changepoint at firsttime point (i.e., the gene responds immediately) and a convex('cv') shaped trend.

Let’s look for a gene that has a more complex trend. Unconstrainedshapes incpam are denoted by'tp'.¹ We can filter theresults for genes with unconstrained shapes and plot one of them.

results(cpo)%>%filter(shape=="tp")

We plot the first gene in the list.

plot_cpam(cpo,gene_id ="g210")

This selection of'tp' suggests that the trend for thisgene does not conform to one of the simpler shape types thatcpam uses. We can exclude'tp' as an optionand forcecpam to choose among the simpler forms by settingshape_type = "shape2" in theplot_cpamfunction ("shape1", the default, allows the'tp'). This choice can be useful for analysis such asclustering. For example:

plot_cpam(cpo,gene_id ="g210",shape_type ="shape2")

Here a monotonic increasing concave shape (‘micv’) is chosen, and we cansee this trend deviates from the data substantially more that theunconstrained shape. See themanuscript for moredetails on the shape types available incpam.

Next we look for a gene with a changepoint by filtering for geneswith changepoints at the third time point.

results(cpo)%>%filter(cp==3)

Again, we plot the first gene in the list.

plot_cpam(cpo,gene_id ="g013")

Clusters

The results function can be used to generate clusters according toselected filters. Theplot_cluster function can then beused to visualise the clusters. With such a small simulated data set, wedon’t have many genes in each cluster, but we can try a few differentclustering options to get an idea of how the function works.

res<-results(cpo)plot_cluster(cpo, res,changepoints =1,shapes =c("cv"))

There are 19 genes with a concave shape and a changepoint at the firsttime point.

More than one shape or changepoint can be provided. For example:

plot_cluster(cpo, res,changepoints =2,shapes =c("dlin","mdcx"))

There are just four genes with decreasing linear or monotonic decreasingconvex shapes which have a changepoint at the second time point.

Clustering can be further refined based on, for example, the rate atwhich the above transcripts attain their maximum values. We illustrateadvanced refinements such as thiscasestudy.

Many more options

This is just a simple example to get you started. The package hasmany more features and options. Check out the two case studies to seehowcpam can be used to analyse real-world data:

• ArabidopsisCase Study
• HumanEmbryo Case Study

¹tp stands forthinplate which is thetype of spline used for the ‘unconstrained’ curves as defined in themgcv package. The curves are still penalised to be smooth,but the shape type is not fixed.