Introduction

This vignette showcases the visualization capabilities of tradeSeq for trajectory-based single-cell analysis. We demonstrate various plotting functions and customization options.

Setup

library(tradeSeq)
library(ggplot2)
library(RColorBrewer)
library(SingleCellExperiment)
library(slingshot)
library(viridis)

# Set seed for reproducibility
set.seed(42)

# Load example data
data(countMatrix, package = "tradeSeq")
data(crv, package = "tradeSeq")
data(celltype, package = "tradeSeq")

counts <- as.matrix(countMatrix)

Fit GAM Models

# Fit GAM models
sce <- fitGAM(counts = counts, sds = crv, nknots = 6, verbose = FALSE)

# Run tests to identify interesting genes
assocRes <- associationTest(sce)
endRes <- diffEndTest(sce)

Gene Expression Smoothers

# Select a significant gene
topGene <- rownames(assocRes[order(assocRes$pvalue), ])[1]

# Basic smoother plot
plotSmoothers(sce, counts, gene = topGene)

# Customized plot with borders and adjusted aesthetics
p <- plotSmoothers(sce, counts, gene = topGene,
                   lwd = 2.5,          # Line width
                   size = 0.8,         # Point size  
                   alpha = 0.6,        # Point transparency
                   border = TRUE) +    # White border around curves
  ggtitle(paste("Expression of", topGene)) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    legend.position = "bottom"
  )
print(p)

# Define custom colors for lineages
customCols <- c("#E74C3C", "#3498DB")

plotSmoothers(sce, counts, gene = topGene,
              curvesCols = customCols,
              border = TRUE,
              alpha = 0.5) +
  ggtitle("Custom Color Scheme") +
  theme_minimal()

# Plot only lineage 1
plotSmoothers(sce, counts, gene = topGene,
              lineagesToPlot = 1,
              curvesCols = c("#2ECC71", "transparent"),
              border = FALSE) +
  ggtitle("Lineage 1 Only") +
  theme_minimal()

Gene Count Plots

plotGeneCount(crv, counts, gene = topGene)

# Show knot positions on the trajectory
plotGeneCount(crv, counts, gene = topGene, models = sce) +
  ggtitle(paste(topGene, "with Knot Positions"))

Multiple Gene Comparison

# Get top 4 significant genes
topGenes <- rownames(assocRes[order(assocRes$pvalue), ])[1:4]

# Create a grid of plots
library(gridExtra)

plots <- lapply(topGenes, function(g) {
  plotSmoothers(sce, counts, gene = g, 
                lwd = 1.5, size = 0.5, alpha = 0.4) +
    ggtitle(g) +
    theme(plot.title = element_text(size = 10, hjust = 0.5),
          legend.position = "none")
})

grid.arrange(grobs = plots, ncol = 2)

Heatmap Visualization

# Get predicted expression for top genes
topGenes <- rownames(assocRes[order(assocRes$pvalue), ])[1:20]

# Predict smooth curves
yhat <- predictSmooth(sce, gene = topGenes, nPoints = 100, tidy = FALSE)

# Scale expression
yhatScaled <- t(scale(t(yhat)))

# Create heatmap using base R
# Check if pheatmap is available
if (requireNamespace("pheatmap", quietly = TRUE)) {
  # Color palette
  colors <- colorRampPalette(c("#3B9AB2", "#78B7C5", "#EBCC2A", "#E1AF00", "#F21A00"))(100)
  
  pheatmap::pheatmap(
    yhatScaled,
    cluster_cols = FALSE,
    cluster_rows = TRUE,
    show_colnames = FALSE,
    color = colors,
    main = "Gene Expression Along Trajectory",
    fontsize_row = 8
  )
}

Knot Evaluation Plot

The evaluateK() function helps determine the optimal number of knots:

# Evaluate different numbers of knots
evalRes <- evaluateK(counts = counts, 
                     sds = crv, 
                     k = 3:10, 
                     verbose = FALSE)

This produces a diagnostic plot showing AIC values for different knot numbers, helping you choose the optimal complexity for your smoothers.

Volcano Plot of Results

# Create volcano plot from association test results
assocRes$negLogP <- -log10(assocRes$pvalue)
assocRes$gene <- rownames(assocRes)

# Define significance
assocRes$significant <- assocRes$pvalue < 0.05 & abs(assocRes$meanLogFC) > 0.5

ggplot(assocRes, aes(x = meanLogFC, y = negLogP, color = significant)) +
  geom_point(alpha = 0.6, size = 2) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray50") +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("gray50", "#E74C3C")) +
  labs(
    title = "Volcano Plot: Association Test Results",
    x = "Mean Log Fold Change",
    y = "-log10(p-value)",
    color = "Significant"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

Ranked Gene Plot

# Rank genes by significance
assocRes <- assocRes[order(assocRes$pvalue), ]
assocRes$rank <- 1:nrow(assocRes)

# Plot top 30 genes
top30 <- assocRes[1:30, ]

ggplot(top30, aes(x = reorder(gene, -negLogP), y = negLogP)) +
  geom_bar(stat = "identity", fill = "#3498DB", alpha = 0.8) +
  coord_flip() +
  labs(
    title = "Top 30 Genes by Significance",
    x = "Gene",
    y = "-log10(p-value)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.text.y = element_text(size = 8)
  )

Publication-Ready Figures

library(patchwork)

# Create individual plots
p1 <- plotSmoothers(sce, counts, gene = topGenes[1], 
                    border = TRUE, alpha = 0.5) +
  ggtitle(topGenes[1]) +
  theme(legend.position = "none")

p2 <- plotSmoothers(sce, counts, gene = topGenes[2], 
                    border = TRUE, alpha = 0.5) +
  ggtitle(topGenes[2]) +
  theme(legend.position = "none")

p3 <- ggplot(assocRes[1:20, ], aes(x = reorder(gene, -negLogP), y = negLogP)) +
  geom_bar(stat = "identity", fill = "#E74C3C", alpha = 0.8) +
  coord_flip() +
  labs(x = "", y = "-log10(p-value)", title = "Top 20 Genes") +
  theme_minimal()

p4 <- ggplot(assocRes, aes(x = meanLogFC, y = negLogP, color = significant)) +
  geom_point(alpha = 0.5, size = 1.5) +
  scale_color_manual(values = c("gray50", "#3498DB")) +
  labs(x = "Mean Log FC", y = "-log10(p-value)", title = "Volcano Plot") +
  theme_minimal() +
  theme(legend.position = "none")

# Combine plots
(p1 + p2) / (p3 + p4) +
  plot_annotation(
    title = "tradeSeq Analysis Summary",
    theme = theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5))
  )

Tips for Effective Visualization

1. Choose appropriate colors: Use colorblind-friendly palettes
2. Add context: Include titles, labels, and legends
3. Consider density: For large datasets, use transparency or subsampling
4. Export high resolution: Use ggsave() with appropriate DPI for publications

# Save publication-quality figure
ggsave("tradeSeq_figure.pdf", width = 10, height = 8, dpi = 300)
ggsave("tradeSeq_figure.png", width = 10, height = 8, dpi = 300)

Session Info

sessionInfo()

## R version 4.6.0 (2026-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 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.26.so;  LAPACK version 3.12.0
## 
## Random number generation:
##  RNG:     Mersenne-Twister 
##  Normal:  Inversion 
##  Sample:  Rounding 
##  
## 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] patchwork_1.3.2             gridExtra_2.3              
##  [3] viridis_0.6.5               viridisLite_0.4.3          
##  [5] clusterExperiment_2.33.0    ggplot2_4.0.3              
##  [7] mgcv_1.9-4                  nlme_3.1-169               
##  [9] slingshot_2.21.0            TrajectoryUtils_1.21.0     
## [11] princurve_2.1.6             SingleCellExperiment_1.35.1
## [13] SummarizedExperiment_1.43.0 Biobase_2.73.1             
## [15] GenomicRanges_1.65.0        Seqinfo_1.3.0              
## [17] IRanges_2.47.1              S4Vectors_0.51.2           
## [19] BiocGenerics_0.59.3         generics_0.1.4             
## [21] MatrixGenerics_1.25.0       matrixStats_1.5.0          
## [23] RColorBrewer_1.1-3          tradeSeq_1.13.13           
## [25] knitr_1.51                  rmarkdown_2.31             
## 
## loaded via a namespace (and not attached):
##   [1] sys_3.4.3            jsonlite_2.0.0       magrittr_2.0.5      
##   [4] farver_2.1.2         vctrs_0.7.3          locfdr_1.1-8        
##   [7] memoise_2.0.1        htmltools_0.5.9      S4Arrays_1.13.0     
##  [10] progress_1.2.3       Rhdf5lib_2.1.0       rhdf5_2.57.0        
##  [13] SparseArray_1.13.2   sass_0.4.10          bslib_0.11.0        
##  [16] plyr_1.8.9           cachem_1.1.0         uuid_1.2-2          
##  [19] buildtools_1.0.0     igraph_2.3.1         lifecycle_1.0.5     
##  [22] iterators_1.0.14     pkgconfig_2.0.3      rsvd_1.0.5          
##  [25] Matrix_1.7-5         R6_2.6.1             fastmap_1.2.0       
##  [28] digest_0.6.39        colorspace_2.1-2     AnnotationDbi_1.75.0
##  [31] phylobase_0.8.12     irlba_2.3.7          RSQLite_3.53.1      
##  [34] beachmat_2.29.0      labeling_0.4.3       httr_1.4.8          
##  [37] abind_1.4-8          compiler_4.6.0       rngtools_1.5.2      
##  [40] bit64_4.8.2          withr_3.0.2          doParallel_1.0.17   
##  [43] S7_0.2.2             BiocParallel_1.47.0  DBI_1.3.0           
##  [46] HDF5Array_1.41.0     MASS_7.3-65          DelayedArray_0.39.2 
##  [49] zinbwave_1.35.0      tools_4.6.0          rncl_0.8.9          
##  [52] ape_5.8-1            glue_1.8.1           h5mread_1.5.0       
##  [55] rhdf5filters_1.25.0  grid_4.6.0           gridBase_0.4-7      
##  [58] cluster_2.1.8.2      reshape2_1.4.5       ade4_1.7-24         
##  [61] gtable_0.3.6         tidyr_1.3.2          hms_1.1.4           
##  [64] BiocSingular_1.29.0  ScaledMatrix_1.21.0  xml2_1.5.2          
##  [67] XVector_0.53.0       foreach_1.5.2        pillar_1.11.1       
##  [70] stringr_1.6.0        limma_3.69.1         genefilter_1.95.0   
##  [73] softImpute_1.4-3     splines_4.6.0        dplyr_1.2.1         
##  [76] lattice_0.22-9       survival_3.8-6       bit_4.6.0           
##  [79] annotate_1.91.0      tidyselect_1.2.1     registry_0.5-1      
##  [82] locfit_1.5-9.12      Biostrings_2.81.2    maketools_1.3.2     
##  [85] pbapply_1.7-4        edgeR_4.11.1         xfun_0.57           
##  [88] statmod_1.5.2        NMF_0.28             pheatmap_1.0.13     
##  [91] stringi_1.8.7        yaml_2.3.12          evaluate_1.0.5      
##  [94] codetools_0.2-20     kernlab_0.9-33       tibble_3.3.1        
##  [97] BiocManager_1.30.27  cli_3.6.6            xtable_1.8-8        
## [100] jquerylib_0.1.4      Rcpp_1.1.1-1.1       png_0.1-9           
## [103] XML_3.99-0.23        parallel_4.6.0       RNeXML_2.4.11       
## [106] blob_1.3.0           prettyunits_1.2.0    scales_1.4.0        
## [109] purrr_1.2.2          crayon_1.5.3         rlang_1.2.0         
## [112] KEGGREST_1.53.0