bio-differential-expression-timeseries-de

name	bio-differential-expression-timeseries-de
description	Analyze time-series RNA-seq data using limma voom with splines, maSigPro, and ImpulseDE2. Identify genes with dynamic expression patterns.
tool_type	r
primary_tool	limma, maSigPro

Time-Series Differential Expression

Identify genes with significant temporal expression patterns in time-course experiments.

Approaches

Method	Best For
limma with splines	Smooth temporal patterns
maSigPro	Multiple time points, regression
ImpulseDE2	Impulse-like patterns
DESeq2 LRT	Discrete time comparisons

limma with Splines

Setup

library(limma)
library(edgeR)
library(splines)

# Load count data
counts <- read.table('counts.txt', header=TRUE, row.names=1)
metadata <- read.table('metadata.txt', header=TRUE)

# metadata should have: sample, time, condition, replicate

Basic Time-Series Model

# Create DGEList
dge <- DGEList(counts=counts)
dge <- calcNormFactors(dge)

# Filter low counts
keep <- filterByExpr(dge, group=metadata$condition)
dge <- dge[keep, , keep.lib.sizes=FALSE]

# Design with natural splines
time <- metadata$time
design <- model.matrix(~ ns(time, df=3))

# voom transformation
v <- voom(dge, design, plot=TRUE)

# Fit model
fit <- lmFit(v, design)
fit <- eBayes(fit)

# Test for any temporal effect (all spline terms)
results <- topTable(fit, coef number

Two Conditions Over Time

# Design for condition-specific time effects
condition <- factor(metadata$condition)
time <- metadata$time

# Interaction model
design <- model.matrix(~ condition * ns(time, df=3))

v <- voom(dge, design, plot=TRUE)
fit <- lmFit(v, design)
fit <- eBayes(fit)

# Genes with different temporal patterns between conditions
# Test interaction terms
results_interaction <- topTable(fit, coef=grep(':', colnames(design)), number=Inf)

Contrasts for Specific Comparisons

# Compare time points within condition
design <- model.matrix(~ 0 + condition:factor(time))
colnames(design) <- gsub(':', '_', colnames(design))

v <- voom(dge, design)
fit <- lmFit(v, design)

# Contrast: Treated_T2 vs Treated_T0
contrast <- makeContrasts(
    early_response = ConditionTreated_time2 - ConditionTreated_time0,
    late_response = ConditionTreated_time6 - ConditionTreated_time0,
    levels = design
)

fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)
results <- topTable(fit2, coef='early_response', number

maSigPro

Specialized for time-series analysis.

Installation

BiocManager::install('maSigPro')

Two-Step Regression

library(maSigPro)

# Create experimental design
# Time, Replicate, Group columns required
edesign <- data.frame(
    Time = metadata$time,
    Replicate = metadata$replicate,
    Control = as.numeric(metadata$condition == 'Control'),
    Treatment = as.numeric(metadata$condition == 'Treatment')
)
rownames(edesign) <- metadata$sample

# Normalize counts
dge <- DGEList(counts=counts)
dge <- calcNormFactors(dge)
norm_counts <- cpm(dge, log=TRUE)

# Create design matrix for polynomial regression
design <- make.design.matrix(edesign, degree=3


fit  p.vectornorm_counts design Q MT.adjust


tstep  T.fitfit step.method alfa


sigs  get.siggeneststep rsq vars


see.genessigssig.genes show.fit disdesigndis
          cluster.method k

Cluster Visualization

# Plot specific clusters
pdf('timeseries_clusters.pdf', width=12, height=10)
see.genes(sigs$sig.genes, show.fit=TRUE, dis=design$dis,
          cluster.method='hclust', k=9,
          newX11=FALSE)
dev.off()

# Get genes per cluster
cluster_genes <- sigs$sig.genes$sig.profiles

ImpulseDE2

For impulse-like expression patterns.

Installation

BiocManager::install('ImpulseDE2')

Run ImpulseDE2

library(ImpulseDE2)
library(DESeq2)

# Create annotation
dfAnnotation <- data.frame(
    Sample = colnames(counts),
    Time = metadata$time,
    Condition = metadata$condition,
    Batch = metadata$batch
)

# Run ImpulseDE2
impulse_results <- runImpulseDE2(
    matCountData = as.matrix(counts),
    dfAnnotation = dfAnnotation,
    boolCaseCtrl = TRUE,
    vecConfounders = c('Batch'),
    scaNProc = 4
)

# Get significant genes
sig_genes <- impulse_results$dfImpulseDE2Results[
    impulse_results$dfImpulseDE2Results$padj < 0.05, ]

DESeq2 Likelihood Ratio Test

For discrete time points without assuming smooth curves.

library(DESeq2)

# Design with time as factor
dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData = metadata,
    design = ~ condition + time + condition:time
)

# LRT: test if time has any effect
dds_lrt <- DESeq(dds, test='LRT', reduced = ~ condition)
results_lrt <- results(dds_lrt)

# Genes with significant temporal pattern
sig_time <- results_lrt[results_lrt$padj < 0.05 & !is.na(results_lrt$padj), ]

Visualization

Expression Profiles

library(ggplot2)

plot_gene_timeseries <- function(gene, counts, metadata) {
    gene_data <- data.frame(
        time = metadata$time,
        condition = metadata$condition,
        expression = as.numeric(counts[gene, ])
    )

    ggplot(gene_data, aes(time, expression, color = condition, group = condition)) +
        geom_point(size = 2) +
        geom_smooth(method = 'loess', se = TRUE, alpha = 0.2) +
        labs(title  gene x   y   
        theme_bw



top_genes  headrownamesresultsorderresultsadj.P.Val 
plots  lapplytop_genes plot_gene_timeseries counts  norm_counts metadata  metadata
librarypatchwork
wrap_plotsplots ncol  

Heatmap with Time Order

library(pheatmap)

# Get significant genes
sig_genes <- rownames(results)[results$adj.P.Val < 0.05]

# Order samples by time
sample_order <- order(metadata$time, metadata$condition)
mat <- norm_counts[sig_genes, sample_order]

# Scale rows
mat_scaled <- t(scale(t(mat)))

# Annotation
anno_col <- data.frame(
    Time = metadata$time[sample_order],
    Condition = metadata$condition[sample_order],
    row.names = colnames(mat)
)

pheatmap(mat_scaled,
         annotation_col = anno_col,
         cluster_cols = FALSE
         show_rownames  
         color  colorRampPalette  

Complete Workflow

library(limma)
library(edgeR)
library(splines)
library(maSigPro)

# Load data
counts <- read.table('counts.txt', header=TRUE, row.names=1)
metadata <- read.table('metadata.txt', header=TRUE, row.names=1)

# Normalize
dge <- DGEList(counts=counts)
dge <- calcNormFactors(dge)
keep <- filterByExpr(dge, group=metadata$condition)
dge <- dge[keep, , keep.lib.sizes=FALSE]
norm_counts <- cpm(dge, log=TRUE


design  model.matrix metadatacondition  nsmetadatatime df
v  voomdge design plot
fit  lmFitv design
fit  eBayesfit


interaction_terms  grep colnamesdesign
results_limma  topTablefit coefinteraction_terms number
sig_limma  rownamesresults_limmaresults_limmaadj.P.Val  


edesign  data.frame
    Time  metadatatime
    Replicate  nrowmetadata
    Control  metadatacondition  
    Treatment  metadatacondition  

rownamesedesign  rownamesmetadata

design_masig  make.design.matrixedesign degree
fit_masig  p.vectornorm_counts design_masig Q
tstep  T.fitfit_masig step.method
sigs  get.siggeneststep rsq vars


all_sig  unionsig_limma rownamessigssig.genessig.profiles
cat all_sig 

Related Skills

• differential-expression/deseq2-basics - Standard DE analysis
• differential-expression/de-visualization - Visualize results
• differential-expression/batch-correction - Handle batch effects
• pathway-analysis/go-enrichment - Functional analysis of clusters