---
title: "Visualization Gallery"
author: "Zaoqu Liu"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: flatly
    highlight: tango
    code_folding: show
vignette: >
  %\VignetteIndexEntry{Visualization Gallery}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 9,
  fig.height = 6,
  fig.align = "center",
  dpi = 150
)
```

# Introduction

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

# Setup

```{r load-packages}
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

```{r fit-gam}
# 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

## Basic Smoother Plot

The `plotSmoothers()` function visualizes the fitted expression curves for a gene.
```{r basic-smoother, fig.height=5}
# Select a significant gene
topGene <- rownames(assocRes[order(assocRes$pvalue), ])[1]

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

## Customized Smoother Plot

```{r custom-smoother, fig.height=5}
# 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)
```

## Custom Color Schemes

```{r custom-colors, fig.height=5}
# 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 Specific Lineages

```{r specific-lineage, fig.height=5}
# 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

## Basic Gene Count Plot

The `plotGeneCount()` function shows gene expression overlaid on the trajectory.

```{r basic-genecount, fig.height=5}
plotGeneCount(crv, counts, gene = topGene)
```

## With Model Knots

```{r genecount-knots, fig.height=5}
# Show knot positions on the trajectory
plotGeneCount(crv, counts, gene = topGene, models = sce) +
  ggtitle(paste(topGene, "with Knot Positions"))
```

# Multiple Gene Comparison

## Comparing Expression Patterns

```{r multi-gene, fig.height=8, fig.width=10}
# 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

## Expression Heatmap Along Trajectory

```{r heatmap-prep}
# 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)
```

```{r heatmap-plot, fig.height=8, fig.width=10}
# 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:

```{r eval-k, eval=FALSE}
# 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

```{r volcano, fig.height=6}
# 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

```{r ranked-genes, fig.height=5}
# 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

## Combined Figure

```{r pub-figure, fig.height=10, fig.width=12}
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

```{r save-example, eval=FALSE}
# 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

```{r session-info}
sessionInfo()
```